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ABSTRACT 


This  project  was  a  three  year  team  effort  (1997  -  2000),  led  by  CFD  Research  Corporation 
(CFDRC),  with  Carnegie  Mellon  University  (CMU),  Georgia  Institute  of  Technology  (GT),  and 
University  of  Florida  (UF).  As  a  result  of  this  project,  we  have  developed,  demonstrated,  and 
validated  a  Software  System  for  Automated  Generation  of  Reduced  or  Compact  Models  of 
Microdevices  from  High  Fidelity  3-D  Simulations,  including:  •  new  concepts  of 
reduced/compact  models  of  electro-mechanical,  fluidic,  and  thermal  phenomena  in  Micro- 
Electro-Mechanical  Systems  (MEMS);  •  interfaces  and  import/export  procedures  for  CAD/EDA 
standards,  to  enable  easy  and  automatic  building  of  three-dimensional  (3-D)  models  from  EDA 
layouts  and  process  data;  •  new,  updated,  and  validated  features  in  CFD-ACE+  for  reliable  high- 
fidelity  simulations  of  coupled  physical  phenomena  specific  for  microdevices;  •  procedures  and 
user-friendly  interfaces  for  extraction  of  compact-model  parameters  (characteristics)  from  3-D 
simulations;  •  output  formats  compatible  with  system-level  simulators  like  SPICE, 
SABER/MAST,  VHDL-AMS;  •  libraries  of  reduced  models  for  MEMS  CAD  system-level 
tools,  like  NODAS  from  CMU;  •  tools  in  CFD-ACE+  allowing  for  use  of  reduced  or  mixed 
dimensionality  in  ID,  2D,  or  3-D  simulations  of  MEMS  devices  and  systems. 

The  most  important  accomplishments  include: 

-  Compact  models  (SPICE  and  SABER)  of  air-damping  in  MEMS,  validated  with  CFD-ACE+ 
3-D  simulations,  and  inserted  into  the  NODAS  MEMS-design-system  library  from  CMU. 

-  Automatic  generation  of  reduced  thermal  models  (R-networks)  for  MEMS  micro-heaters. 

-  Developed  mixed-dimensional  and  compact  fluidic  models  of  microchannels,  microvalves, 
micropumps,  droplet  generators,  and  synthetic  jets,  and  procedures  of  their  generation. 

-  A  new  software,  CFD-Micromesh,  for  automatic  3-D  solid  model  generation  and  meshing 
from  MEMS  layouts  imported  directly  from  CAD/EDA  systems.  The  popular  layout  formats 
(GDSII,  CIF,  DXF)  are  imported  easily,  allowing  direct  coupling  with  several  commercial  IC 
design  tools. 

-  At  CMU:  *  The  compact  squeeze  film  model  was  compared  with  experimental  data  from  a 
CMOS-MEMS  bandpass  filter.  The  NODAS  simulation  results  with  new  model  are  in 
excellent  agreement  with  the  experimental  measurements  (error  decreased  from  20%  to  2%). 
*  A  first-order  model  for  lateral  damping  within  NODAS.  *  A  new  data  format  GBV 
(Geometry,  Boundary  and  Volume  conditions)  for  interchange  of  simulation  data  between 
Cadence  and  CFDRC  simulators. 

-  At  GT:  *  Fabricated  synthetic  jets  with  integrated  MEMS  modulators  and  measured  with 
PIV.  *  Measured  characteristics  of  arrays  of  synthetic  jets  in  active  control  of  microsystem 
temperature.  *  The  data  measured  with  infra-red  camera  was  successfully  compared  with 
numerical  simulations  at  CFDRC,  including  the  use  of  reduced  models  in  jet  arrays. 

-  At  UF:  *  New  matrix  reduction  techniques  for  structure  dynamics,  based  on  Lanczos 
method,  WYD  Ritz  algorithm.  *  Validation  of  the  implementation  of  Lanczos  algorithm 
through  comparison  with  a  MEMS-bridge  structure  from  MIT.  *  Testing  the  new  methods  on 
the  crab-leg  flexure  MEMS  structure.  *  The  tests  showed  from  600x  to  10,000x  acceleration 
in  comparison  to  full  3-D  FEM  computation  time.  *  Implementation  of  the  WYD  algorithm  in 
the  CFD-ACE+  code  of  CFDRC. 

The  high-fidelity  and  compact  models  developed  in  this  effort  have  already  been  used  in  design 
work  by  Kodak  (ink-jets);  Delphi  Automotive  (air-damping  in  accelerometers).  Lucent 


(micromirrors),  and  contributed  to  the  DARPA  HERETIC  Program  in  modeling  and  design  of 
active  cooling  of  VCSEL  Arrays  {Georgia  Tech,  Honeywell),  Microprocessor  Chips  {Georgia 
Tech,  Intel).  This  effort  resulted  also  in  awarding  CFDRC  with  a  new  SBIR  contract  “High 
Speed  Inertial  MEMS  Sensors  with  Field  Emitter  Array  Readout”  from  U.S.  Air  Force 
(AFRLA^SSE)  to  design  novel  accelerometers  and  gyroscopes  for  satellite  and  vehicle 
positioning. 
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INTRODUCTION 


1.1  Background 

This  project  began  in  August  1997.  It  was  a  three  year  team  effort  between  CFD  Research 
Corporation,  Carnegie  Mellon  University,  Georgia  Institute  of  Technology,  and  University  of 
Florida. 

The  objective  of  this  project  was  to  develop,  demonstrate  and  validate  a  software  system  for 
automated  generation  and  selection  of  reduced  compact  models  (CMs)  of  MEMS  devices  from 
high  fidelity  3-D  models.  The  developed,  extracted,  and  validated  CMs  would  be  used  in  two 
types  of  environments: 

•  within  a  system-level  behavior  simulation  (like  SPICE  or  SABER)  to  support  mixed 
technology  CAD  tools  for  microsystems;  and 

•  as  subgrid  models  within  CFD-ACE+,  a  3-D  field  simulation  environment  used  for  larger 
scale  (system/package)  analysis. 

Major  activities  in  this  program  were  organized  to  take  advantage  of  the  expertise  of  each  team 
member: 

•  CFDRC:  (Dr.  Marek  Turowski,  Dr.  Andrzej  Przekwas,  Dr.  Zhijian  Chen)  project 
coordination,  new  reduced/compact  model  concepts,  software  development,  validation, 
demonstration,  and  applications. 

•  Carnegie  Mellon  University:  (Dr.  Gary  Redder,  Dr.  Tamal  Mukherjee)  reduced  model 
generation  for  inertial  sensors,  interfaces  for  geometry  extraction  for  3-D  modeling  from 
MEMS  standard  CAD/EDA  systems,  interfaces  to  MEMS  synthesis  database,  library  of 
reduced  models,  software  validation. 

•  Georgia  Tech:  (Dr.  Mark  Allen,  Dr.  Ari  Glezer)  fabrication  of  batch  of  synthetic 
microjet  devices,  full  field  experimental  characterization  of  microjet  devices, 
experimental  extraction  of  compact  models  and  verification  against  computationally 
generated  models,  development  of  fabrication  dependent  multi-disciplinary  material 
property  database. 

•  University  of  Florida:  (Dr.  Loc  Vu-Quoc)  development  of  formal  mathematical  methods 
for  reduced  model  generation  from  FEM  simulators,  development  and  implementation  of 
model  extraction  and  reduction  capability  in  FEM-STRESS  and  CFD-ACE-t-  for 
stress/deformation  simulations  in  MEMS  high-level  design. 

This  report  documents  the  work  performed  during  the  three  years  of  the  project. 
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1.2  Project  Motivation 


Micro-Electro-Mechanical  Systems  (MEMS)  consist  of  interconnected  arrays  of  mechanical, 
electrostatic,  thermal,  structural,  optical,  fluidic,  magnetic,  and  electronic  components.  Individual 
components  can  be  simulated  using  multi-dimensional  high-fidelity  computational  tools  based  on 
first  principles  and  Partial  Differential  Equations  (PDEs).  For  system  level  design,  where  large 
numbers  of  devices  are  interconnected,  use  of  detailed  3-D  models  of  each  component  is 
unrealistic.  A  more  practical  approach  is  to  develop  simplified  but  accurate  enough  submodels. 
Reduced  or  Compact  Models  (CMs)  expressed  in  terms  of  Ordinary  Differential  Equations 
(ODEs)  or  ID-PDEs  would  perform  this  task.  For  example,  a  3-D  Navier-Stokes  equation 
system  could  be  reduced  to  ID  Stokes  or  Couette  type  parametric  relation  CMs,  mixing  reaction 
channels/chambers  could  be  reduced  to  ID  plug  flow  or  OD  stirred  reactor  equivalents. 

Neither  Electronic  Design  Automation  (EDA)  nor  mechanical  CAD  can  adequately  support 
chips  that  incorporate  MEMS  devices.  The  scale  is  much  different  and  the  materials  properties 
are  different.  At  the  beginning  of  the  DARPA  Composite  CAD  Program  there  were  many 
challenges  facing  those  developing  design  tools  for  MEMS-on-a-chip.  These  included  a  need  to 
integrate  EDA  and  TCAD  (Technology  CAD)  tools  with  solvers  for  mechanical,  electrostatic, 
magnetic,  thermal,  fluidic  and  optical  domains;  the  slow  speed  of  3-D  mechanical  finite-element 
analysis  tools;  a  lack  of  software  models;  and  a  need  to  do  extensive  physical  prototyping.  The 
development  of  design-synthesis  and  design-validation  tools  would  greatly  simplify  the  process. 

Similarly  to  the  VLSI  electronic  industry,  there  is  a  great  need  of  "MEMS  synthesis"  software 
tools.  A  MEMS  synthesis  tool  would  generate  a  physical  layout  cell  for  an  on-chip  MEMS 
device  from  user-defined  parameters.  Synthesis  process  needs  fast  simulation  tools  to  optimize 
the  design.  This  requires  reduced  or  compact  models  of  MEMS  devices,  based  on  algebraic 
models,  not  partial  differential  equations.  In  such  a  case,  the  speed  of  simulation  is  exchanged  for 
accuracy  of  modeling.  There  was  a  big  need  to  make  adequate  reduced  models  to  ensure  that 
design  constraints  have  been  met.  Model  generation  has  been  a  key  problem  area  in  MEMS 
synthesis  and  system-level  design. 

There  was  also  a  great  need  to  develop  design  aids  for  a  "micro-bio-chemical  network”.  In  brief, 
such  a  system  would  include  "wells"  that  could  contain  reaction  or  observation  chambers,  and 
"channels"  that  might  transport  microliters  of  fluids  between  wells.  Researchers  are  hoping  to  put 
hundreds  of  wells  and  channels  on  a  single  chip,  making  it  possible  to  analyze  a  number  of 
different  reagents  from  very  small  samples.  To  be  able  to  design  an  entire  system  like  this,  there 
was  a  need  of  software  that  can  use  the  results  from  a  computational  fluid  dynamics  (CFD) 
solver,  and  allow  designers  to  run  simulations  similar  in  nature  to  SPICE.  The  goal  is  a  tool  that 
can  generate  and  lay  out  microfluidic  circuits  according  to  the  user's  constraints.  Another  task 
was  to  take  the  results  of  a  full  3-D  simulation,  generate  a  behavioral  model  and  bring  it  into 
SPICE  for  an  analysis  of  interactions  with  electrical  systems. 

The  existing  EDA  software  needed  to  be  extended  in  three  areas:  physical  design,  analysis  and 
behavioral  design.  In  the  physical-design  area,  there  was  a  great  need  of  a  solid  model  generator 
using  directly  layouts  from  EDA  layout  editors,  which  would  produce  a  3-D  representation  of 
what  a  mixed-technology  chip  will  look  like.  In  the  analysis  phase,  that  representation  has  to  be 
passed  to  a  mixed-domain  boundary-element  or  finite-element  analysis  tool.  This  tool  should  be 
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able  to  co-simulate  in  the  electrostatic,  thermal,  and  fluidic  domains.  For  the  behavioral  level,  a 
behavioral  model  builder  had  to  be  developed  so  that  SPICE  type  simulations  tools  should  be 
able  to  handle  mechanical,  thermal,  and  fluidic  variables. 

1.3  Project  Objectives 

The  objective  of  this  project  was  to  develop,  demonstrate,  and  validate: 

•  Software  System  (in  CFD-ACE+)  for  Automated  Generation  of  Reduced  or 
Compact  Models  of  Microdevices  from  High  Fidelity  3-D  Simulations 

including: 

new  concepts  of  reduced/compact  models  of  electro-mechanical,  fluidic,  and  thermal 
phenomena  in  MEMS; 

interfaces  and  import/export  procedures  for  CAD/EDA  standards,  to  enable  easy  and 
automatic  building  of  3-D  models  from  EDA  layouts  and  process  data; 

-  new,  updated,  and  validated  features  in  CFD-ACE+  for  reliable  high-fidelity 
simulations  of  eoupled  physical  phenomena  specific  for  microdevices; 

-  procedures  and  user-friendly  interfaces  for  extraction  of  compact-model  parameters 
(characteristics)  from  3-D  simulations; 

-  output  formats  compatible  with  system-level  models,  for  simulators  like  SPICE, 
SABER/MAST,  VHDL-AMS; 

-  libraries  of  reduced  models  for  MEMS  CAD  system-level  tools,  like  NODAS; 

-  tools  in  CFD-ACE-i-  allowing  for  use  of  reduced  or  mixed  dimensionality  in  ID,  2D, 
or  3-D  simulations  of  MEMS  devices  and  systems. 

Classes  of  Reduced  Models  to  Be  Considered: 

A.  Phenomenological  (behavioral),  based  on  fitted  curves,  lookup  tables,  etc., 

-  simulators:  SABER  (MAST,  VHDL-AMS),  SPICE,  Spectre  (from  Cadence); 

B.  Parametric  Models  based  on  Physical  Decomposition, 

-  physics-based  decomposition  of  the  problem  into  sub-problems  that  can  be  described 
by  analytical,  parameterized  equations, 

-  simulators:  SABER,  HSPICE; 

C.  Discrete  -  reduced-size  physical  distributed  models  based  on  PDEs, 

-  "mixed  mode"  simulators:  2D/3-D  PDE  solvers  coupled  with  external  system/circuit; 

D.  Formally  Reduced  Matrix  Systems, 

-  FEM  models  reduced  with  Wilson-Ritz,  Lanczos,  or  Amoldi  approach; 

-  simulators:  matrix/vector  solvers  coupled  with  external  system/circuit  (e.g. 
MATLAB). 

Domains  and  Devices  Identified  for  Investigation  in  This  Project: 

1.  Fluidic 
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-  damping:  inertial  sensors,  plates,  plates  with  holes,  combs,  cantilever 
beams/plates,  moving  micromirrors; 

-  microjets:  single  jet,  multiple  jets,  multiple  jets  for  electronic  cooling,  multiple  jets 
for  flow  control; 

-  microchannels 

2.  Electro-Mechanical 

-  beams,  plates,  bridges  actuated  electrostatically 

-  electrostatic  inertial  sensors  (comb  drives  and  resonators) 

3.  Thermal 

-  microheaters 

-  MEMS  packages 

-  synthetic-jet  based  cooling  systems 

1.4  Summary  of  Accomplishments 

All  of  the  accomplishments  aimed  for  fulfillment  of  the  main  objective  of  this  project  and  its 
components,  as  described  in  the  previous  section.  Following  subsections  list  the  most  important 
accomplishments  categorized  by  their  role  in  the  contribution  to  the  main  objective. 


New  Reduced/Compact  Models  of  Microdevices 

•  New  nonlinear  compact  models  of  squeeze  film  damping  in  MEMS,  validated  with  CFD- 
ACE+  3-D  simulations,  and  implemented  in  SPICE  and  SABER  fonnats  (Figure  1),  and 
implemented  in  NODAS  MEMS-CAD  system  at  CMU. 


SPICE  Model  -  Equivalent  Circuit 
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Figure  1.  Compact  models  of  squeeze  film  damping  for  system-level  simulators  (SPICE  and 
SABER),  derived  from  high-fidelity  simulations. 
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•  The  new  compact  squeeze  film  model  was  successfully  compared  with  experimental  data 
from  a  CMOS-MEMS  bandpass  filter.  The  NOD  AS  simulation  results  with  the  new 
damping  model  are  in  excellent  agreement  with  the  measurements. 

•  New  compact  model  of  lateral  (shear)  damping  in  MEMS,  validated  with  CFD-ACE+ 
(Figure  2). 


Figure  2.  Equivalent  circuit  representing  the  lateral  shear  damping  forces,  derived  and  validated 

with  CFD-ACE+  high-fidelity  simulations. 


•  Developed  mixed-dimensional  and  reduced  fluidic  models  of  microchannels, 
microvalves,  micropumps,  droplet  generators,  synthetic  jets,  and  piezoelectric 
micropump. 

•  Novel  concept  of  compact  models  for  synthetic  jets,  using  a  polyhedral  control  volume 
capability  of  CFD-ACE-t-  software  tool,  to  model  complex  dynamic  3-D  shapes  with 
moving  walls  with  multiple  inlets/outlets  with  a  single  cell  "super  element".  The  model 
was  validated  against  3-D  high-fidelity  simulation  data  obtained  for  a  range  of 
parameters. 

•  New  equivalent-circuit  models,  in  SPICE  and  SABER  formats,  for  arbitrary- shaped 
fluidic  microchannels,  and  procedures  of  their  generation  (Figure  3). 
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Figure  3.  CFD-ACE+  high-fidelity  simulations  (a  and  b)  are  used  to  generate  compact  models  of  an 
arbitrary  nonlinear  fluidic  microchannel,  for  SPICE  (c)  or  for  SABER  (d). 

•  New  Matrix  Reduction  Techniques  for  structure  dynamics:  Lanczos  method  and  WYD 
Ritz  algorithm,  developed  at  UF.  Tests  showed  up  to  10,000x  acceleration  in  comparison 
to  full  3-D  FEM  computation  time. 


•  Enhancement  of  CFDRC  software  to  process  CAD  data  for  MEMS  design.  We  are  now 
able  to  read  input  files  in  the  formats:  DXF,  IGES,  PATRAN,  CIF,  GDS  11,  and  several 
others. 

•  Two  new  data  formats  have  been  implemented  in  CFD-GEOM,  enabling  input  and 
conversion  from  external  files  coming  from  Electronic  Design  Automation  (EDA)  tools: 

*  CIF  (Caltech  Intermediate  Form) 

*  GDS  II  (Calma  stream  format) 

•  Implementation  of  VRML  (Virtual  Reality  Modeling  Language)  into  the  computational 
environment  of  CFDRC.  An  interface  between  CFDRC's  FastBEM  program  and  VRML- 
based  GUI  have  been  developed. 

•  A  new  software,  CFD-Micromesh,  has  been  developed  during  this  project  at  CFDRC, 
for  automatic  3-D  solid  model  generation  and  meshing  from  MEMS  layouts  imported 
directly  from  CAD/EDA  systems.  The  popular  layout  description  formats,  CIF  and 
GDSII,  are  imported  easily  into  CFD-Micromesh,  allowing  coupling  the  software  with 


several  commercial  IC  design  tools.  The  new  tool  is  fast  and  equipped  with  user-friendly 
graphical  user  interface  (GUI)  -  see  Figure  4. 


Figure  4.  CFD-Micromesh:  fast,  automatic  generation  of  3-D  device  models  from  layouts  (EDA 

tools)  and  process  data. 

By  clicking  only  a  single  button  in  CFD-Micromesh,  a  three-dimensional  solid  model  is 
generated  fully  automatically  from  layout,  as  shown  in  Figure  4.  Building  the  3-D  model  as 
well  as  automatic  generation  of  a  full  3-D  unstructured  computational  mesh,  again  only  by 
clicking  a  single  button  in  CFD-Micromesh,  is  performed  typically  within  several  minutes 
on  a  current  PC  workstation.  The  automatically  generated  three-dimensional  device  model, 
with  the  unstructured  3-D  computational  mesh  in  DTF  format,  is  imported  directly  into 
CFDRC's  simulator  CFD-ACE+. 

Developed  interfaces  between  CFDRC  software/databases  and  the  following  CAD  tools: 

*  SABER 

*  CADENCE 

An  agreement  with  Cadence  Design  Systems,  Inc.  has  been  negotiated  to  develop  interfaces 
and  integration  of  CFDRC's  MEMS  simulation  software  with  Cadence  CAD  tools  (see 
Figure  5).  See  Section  3.1  for  more  details. 
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Figure  5.  CFD-Micromesh  is  launched  directly  from  Virtuoso  Layout  Editor  in  Cadence  Design 

Framework  for  3-D  Modeling  and  Visualization. 


A  new  data  format  (called  GBV  —  Geometry,  Boundary  and  Volume  conditions)  for 
interchange  of  simulation  data  between  Cadence  Design  Framework  H  schematic  view 
and  CFDRC  ACE+  simulator  has  been  developed  at  CMU. 


New  Features  in  CFD-ACE+  for  High-Fidelitv  Simulations  Specific  for  MEMS 

•  Second  order  boundary  conditions  for  wall  treatment. 

•  Second  order  time  treatment,  through  modified  Crank-Nicolson  scheme. 

•  Automatic  grid  remeshing  for  arbitrary  moving  boundaries/bodies  in  6  degrees  of 
freedom  (DOFs).  A  trans-finite  interpolation  (TFI)  technique  with  mesh  smoothing  was 
developed  and  implemented,  which  allows  for  effective  simulation  of  dynamic  3-D 
motion  of  mirrors,  membranes,  gyroscopes,  and  other  MEMS  structures. 
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•  New  capability  in  CFD-ACE+  for  extraction  lumped  parameters  (normal/tangential  force, 
torque,  total  pressure  and/or  flow)  for  compact  model  generation. 

•  "Property  Manager"  -  new  interactive  framework  software  for  MEMS  material  property 
database. 

Procedures  And  User-Friendly  Interfaces  For  Extraction  Of  Compact-Model  Parameters 

fCharacteristics)  From  High-Fidelity  Simulations 

•  Automatic  generation  of  reduced  thermal  models  (R-networks)  for  MEMS  micro-heaters. 

•  Capability  to  perform  parametric  runs  in  CFD-ACE+  through  graphical  user  interface,  for 
mechanical,  electrostatic,  and  fluidic  simulations  both  steady-state  and  transient. 

•  A  dedicated  computer  program  Squeeze,  written  in  portable  C++  language,  for  automatic 
generation  of  input  list  in  SABER  or  SPICE  format,  describing  equivalent-circuit  model 
of  squeeze-film  damping  between  moving  plates. 

•  Procedures  for  automatic  generation  of  compact  model  (extraction  of  lumped 
capacitances)  of  a  comb-drive  resonator,  using  results  of  high-fidelity  3-D  electrostatic 
simulations  with  FastBEM  from  CFDRC. 

•  Procedures  for  automatic  generation  of  equivalent-circuit  models  of  microfluidic  devices 
from  parametric  runs  of  CFD-ACE+;  examples  of  SPICE  and  SABER  compact  models 
of  Tesla  valve. 

•  Procedures  of  FEM  model  reduction  based  on  Amoldi  method  for  electro-mechanical 
microdevices,  developed  at  Prof.  Jacob  White’s  group  at  MIT  within  the  Composite  CAD 
Program,  have  been  implemented  and  installed  at  CFDRC  and  demonstrated  and  tested  in 
coupling  with  CFD-ACE+  code. 

Libraries  Of  Reduced  Models  for  MEMS  CAD  Svstem-Level  Tools 

•  The  new  compact  models  of  squeeze  film  damping  have  been  incorporated  into  the 
NODAS  MEMS-CAD  system  library  from  CMU. 

•  The  new  compact  models  of  lateral  shear  damping  have  also  been  inserted  into  the 
NODAS  MEMS-CAD  system  library. 

•  The  frequency  dependent  transfer  function  of  lateral  damping  has  been  implemented  in 
Verilog-A  as  a  Laplace  transform. 


Validation/Demonstration  of  the  New  CFD-ACE+  Features  and  the  New  Reduced  Models 


•  At  Georgia  Tech  (GT): 

-  Fabricated  micromachined  synthetic  jets  with  integrated  MEMS  modulators. 
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Several  methods  have  been  used  to  test  the  synthetic  jets  and  modulator  arrays:  visual 
deflection  test,  flow  modulation  test,  Schlieren  visualization,  pitot  tube,  x-wire 
anemometer,  and  particle  image  velocimetry  (PIV). 

Characteristics  of  performance  of  arrays  of  synthetic  jets  in  local  active  control  of 
microsystem  temperature  were  measured  with  infra-red  (IR)  camera. 

•  The  PIV  flow  data  from  GT  were  compared  successfully  with  2D  and  3-D  numerical 
simulations  from  CFD-ACE+. 

•  The  temperature-system  data  measured  at  GT  with  IR  camera  were  compared 
successfully  with  numerical  simulations  at  CFDRC,  including  the  use  of  reduced  models 
in  jet  arrays. 

•  At  University  of  Florida  (UF): 

-  Validation  of  the  implementation  of  Lanczos  algorithm  through  comparison  with  a 
MEMS  structure  from  MIT. 

•  At  Carnegie  Mellon  University  (CMU): 

-  Experimentally  verified  the  squeeze  film  model  using  a  three  resonator  CMOS 
micromachined  mechanical  filter  to  show  that  the  improved  squeeze  film  model 
reduced  error  in  quality  factor  (Q)  computation  from  20%  to  2%. 

-  Experimentally  verified  the  lateral  damping  model  using  folded  flexure  MUMPS 
resonators.  Inclusion  of  finite  size  effects  reduced  errors  from  20%  to  8%. 

•  At  CFDRC,  both  the  high-fidelity  3-D  CFD  simulations  of  Tesla  valve  and  its  compact 
models  in  SPICE  and  Saber  were  successfully  compared  with  characteristics  measured  at 
University  of  Washington. 

The  mixed-dimensionality,  multi-energy  domain  capability  in  CFD-ACE-i-  is  the  first  ever 
computational  software  applicable  to  both  micro-scale,  system-scale  and  mixed  micro/system- 
scale  simulations.  The  macro-scale  (package)  allows  the  user  to  verify  the  performance  of  the 
microdevice  (sensor,  transducer,  actuator,  power  device,  optoelectronic  array)  in  an  external 
system  or  environment. 

"Success  Stories": 

The  Project  Team  has  had  a  strong  commitment  to  demonstrating  and  validating  the  usefulness 
of  the  newly  developed  tools  to  real  problems  relevant  to  both  the  defense  industry  and  the 
commercial  market.  Listed  below  are  some  examples  of  such  applications  contributing  to  the 
“Success  Stories”  of  the  DARPA  Composite  CAD  Program: 

•  CFDRC's  High-Fidelity  and  Compact  Models  of  Ink-Jets  used  by  Kodak. 

•  Reduced  Models  of  Synthetic  Jets  used  for  Virtual  Flight  Control  Simulation. 

•  Air-Damping  Compact  Models  used  by  Delphi  Automotive  Systems. 

•  Air-Damping  Compact  Models  under  development  for  Lucent  Technologies,  for  analysis 
and  design  of  MEMS-based  micromirrors  for  fast  optical  switching. 
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•  Reduced  Models  of  Synthetic  Jets  will  be  used  within  the  DARPA  HERETIC  Program 
(Heat  Removal  by  Thermo-Integrated  Circuits),  where  CFDRC  will  support  modeling 
and  design  of  Arrays  of  Synthetic  Microjets  for  active  cooling  of: 

-  VCSEL  Arrays  (Linear  and  2D)  -  Georgia  Tech  +  Honeywell 

-  Microprocessor  Chips  -  Georgia  Tech  +  Intel. 

•  CFDRC  has  been  awarded  a  new  SBIR  contract  in  2000,  for  U.S.  Air  Force 
(AFRUVSSE),  to  develop  “High  Speed  Inertial  MEMS  Sensors  with  Field  Emitter  Array 
Readout”  for  application  in  novel  gyroscopes  for  satellite  and  vehicle  positioning,  etc. 


2  NEW  FEATURES  IN  CFD-ACE-h 

2.1  Automatic  Grid  Remeshing 

A  new  capability  has  been  implemented  in  CFD-ACE-h  at  CFDRC  to  facilitate  automatic  grid 
remeshing  for  arbitrary  moving  boundaries/bodies  in  6  degrees  of  freedom  (DOFs).  A  trans- 
finite  interpolation  (TFI)  technique  with  mesh  smoothing  was  developed.  This  functionality 
allows  effective  simulation  of  dynamic  3-D  motion  of  mirrors,  membranes,  gyroscopes,  and 
other  MEMS  structures.  Details  are  described  below. 

General  Dynamic  Grid  Generation 

A  new  capability  of  the  general  dynamic  grid  generation  using  TFI  (Trans-Finite-Interpolation) 
methodology  has  been  developed  and  implemented  in  CFD-ACE+  code,  to  enhance  the 
numerical  simulation  of  MEMS  flow  problems.  Due  to  the  deformation  of  boundary,  it  is  often 
required  to  remesh  the  grids.  Without  the  technique  of  TFI,  a  user  actually  needs  to  develop 
grid  generation  packages  for  each  individual  geometry.  It  is  sometimes  very  difficult  for  a 
complicated  geometry.  Sometimes  the  transformation  without  TFI  technique  cannot  guarantee 
the  similarity  of  geometry. 

Without  delving  into  details,  TFI  transforms  any  curvilinear  grids  on  ratio  of  arc  lengths  in  each 
direction.  It  is  thus  inherently  “structured”.  Based  on  the  deformations  on  boundaries,  TFI 
creates  a  linear  topology  mapping  to  calculate  deformation  inside  the  domain,  i.e. 


dri  =/((ir^c) 


(2.1) 


where,  dri  is  the  deformation  vector  of  internal  grids  and  dr/,^  is  the  deformation  vector  on  the 
boundary.  TFI  creates  the  relation  function  by  calculating  the  ratio  of  the  arc  length  of  each  grid 
point  to  the  boundary  points  (see  Figure  6). 


f  = 


I 

j=W,E,S,N,L.H 


(2.2) 


where,  Sj  is  the  arc  length  between  point  i  to  j,  and  then  remeshing  interior  grids  by 
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r„  =  To  +  dr^ 


(2.3) 


where,  r„  and  are  new  and  old  meshes.  TFI  preserves  the  similarity  of  geometry  during  the 
transformation  by  linear  mapping  relation  function.  Figure  7  shows  the  results  of  different 
transformations. 


S 

Figure  6.  Grid  point  in  a  structured  mesh  and  the  boundary  points. 


(b)  original  remesh 


Figure  7.  Illustration  of  Dynamic  Grid  Remeshing. 


2.2  Multiple  Parametric  Runs 

To  enable  automatic  generation  of  reduced  models  for  MEMS  applications,  it  is  often  needed 
that  some  parameters  of  several  3-D  simulations  need  to  be  changed  dynamically  for  subsequent 
computations.  The  most  typical  example  is  stepping  through  oscillation  frequency  of  a  device 
simulation,  to  obtain  its  characteristics  in  frequency  domain,  or  stepping  through  a  boundary 
pressure  value  for  a  microfluidic  device  to  obtain  its  static  pressure-flow  characteristic. 

A  new  feature  of  CFD-ACE+,  that  enables  such  multiple  parametric  runs  of  the  program,  was 
developed  and  implemented  in  this  project.  CFD-ACE+  solver  reads  the  definitions  of  changing 
analysis  parameters  from  an  external  text  file  (file  name  with  the  extension  .par)  and 
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automatically  re-runs  the  jobs,  with  new  input  parameters  for  each  run.  The  new  feature  helps  to 
calculate  multiple  cases  without  stopping  the  software.  The  format  of  a  .par  file  looks  as  follows: 

Case  n 

Surface  m  ;  BC_name,  motion_type  =  formula 


where, 

n=  1,2,3..., 
m  =  1, 2,  3  ..., 
BC_name 
motion_type 


-  number  of  the  case; 

-  number  of  the  surface; 

-  name  of  a  specific  boundary; 

-  definition  of  a  motion. 


Below  is  an  example  of  a  .par  file  for  two  automatic  runs  of  simulations,  with  oscillation 
frequencies  100  Hz  and  200  Hz. 


Case  1 

Surface  1  :  side_wall,  Translation_y  =  0.01*sin(2*3.14*100*t) 


Od.  S0  2 

Surface  1  :  side_wall,  Translation_y  =  0 . 01*sin(2*3 . 14*200*t) 


Note  that  a  user  has  a  possibility  to  define  multiple  formulas  for  a  particular  surface,  which  has 
to  be  assigned  a  specific  name  before  the  job  is  run.  The  CFD-ACE+  software  will  read  all  the 
formulas  from  the  .par  file  and  automatically  run  all  the  subsequent  jobs,  saving  required  results 
into  another  disk  file  {.dat).  The  required  results  and  the  output  format  specification  can  be 
defined  in  an  optional  text  file  with  extension  .fint  (see  next  section). 

A  new  graphical  user  interface  panel  has  been  also  introduced  into  CFD-GUI  to  allow  easy  setup 
of  the  parametric  runs  through  graphical  windows.  The  relevant  plots  of  requested  characteristics 
can  be  ordered  through  Simulation  Manager  in  CFD-ACE+. 

2.3  Output  of  Integral  Surface  Results  for  Static  or  Transient  Characteristics 

To  be  able  to  use  directly  the  results  of  CFD-ACE+  simulations  to  plot  device  characteristics,  or 
time-dependent  waveforms  of  selected  variables,  you  can  order  writing  the  selected  integral 
surface  values  to  a  text  file  in  the  current  simulation  directory.  The  file  with  results  (text  tables) 
will  have  the  extension  .dat.  The  data  from  this  file  can  also  be  used  directly  to  create  a 
compact  model  of  the  analyzed  device. 

To  define  which  values  you  want  to  be  printed  in  the  output  table  files,  you  have  to  create  an 
external  text  file,  called  "format  file"  with  extension  .fint.  This  file  should  include  definitions  of 
surface-output  or  integral  quantities,  which  you  want  from  ACE+  simulations,  and  optionally, 
scaling  factors  for  them  (to  obtain  ready  results  in  desired  units). 

Below  is  an  example  of  a  .fint  file  created  for  a  Tesla  valve  simulation,  to  obtain  static  flow- 
pressure  characteristics.  In  the  Tesla.fint  file,  we  order  an  output  text  file  including  the  Mass 
Flow  value  and  the  mean  pressure  value  P_mean  at  the  Channel_Outlet  surface,  and  the  same  for 
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the  Channel_Inlet  surface.  The  Mass  Flow  value  is  scaled  by  (60*2*  le6)  to  obtain  results 
directly  in  |a,L/min  (the  original  value  would  be  kg/s). 


Example  of  format  file:  Tesla.fmt 


★  ★★***•*•★******★*★*★*★**★**■*■★***★***★★**★***★**★*★★***★ 
Surface  :  Channel_Outlet 
Variable  :  Mass  Flow,  M*60*2*le6 
Variable  :  P_mean 


Surface 

Variable 

Variable 


Channel_Inlet 

Mass  Flow,  M*60*2*le6 

P_mean 


****************************************************** 


An  example  of  an  output  file,  for  the  combined  parametric  runs  using  Tesla.par  file,  and  the 
integral  surface  output  data  ordered  with  Tesla.fmt  file,  is  shown  below.  Note  that  the  Mass  Flow 
value  is  scaled,  and  is  printed  directly  in  pL/min,  and  the  P_mean  pressure  values  are  in  pascals 
(N/m^). 


Example  of  file:  Tesla_CHANNEL_OUTLET.dat 


Mass  Flow 
P_mean 


Case 

1 

Iter . 

# 

100 

1.556542E+03 

6.869466E+03 

Case 

2 

Iter. 

# 

100 

2.385028E+03 

1.288921E+04 

Case 

3 

Iter . 

# 

100 

3.044021E+03 

1.857335E+04 

Case 

4 

Iter . 

# 

100 

3 .610177E+03 

2 .405173E+04 

Case 

5 

Iter. 

# 

100 

4.114709E+03 

2 .938847E+04 

Case 

6 

Iter . 

# 

100 

4.574686E+03 

3 .461645E+04 

★******************************************************** 

The  tabularized  results  from  the  output  file  can  be  imported  directly  into  many  plotting 
programs,  like  Excel,  Xmgr,  etc.  An  example  of  the  results  imported  into  Excel,  and  the 
characteristic  plotted  directly  from  the  parametric  CFD-ACE+  simulations,  are  shown  below. 
The  computed  results  are  plotted  together  with  experimental  measurements  of  the  Tesla  valve, 
obtained  and  published  by  R.  Bardell  and  F.  Forster  from  University  of  Washington  (ASME 
1998,  DSC  -  Vol.66). 
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Tesla  45 


•  Experiment  -  Forward  Flow 
O-  CFD-ACE+  Forward  Flow 


Figure  8.  Tesla  Valve  characteristics  plotted  directly  from  the  parametric  CFD-ACE+ 
simulations,  together  with  experimental  results  from  the  University  of  Washington. 

2.4  Material  Properties  Database  Manager 

A  reusable,  open  database  for  MEMS/microdevices  materials  data  for  structural,  thermal,  and 
electromagnetic  properties  has  been  developed; 

-  CFDRC  has  developed  "Property  Manager"  -  prototype  interactive  framework  software  for 

MEMS  material  property  database.  Property  Manager  allows  for  very  flexible  material 
property  description  and  storing  of  data,  including  searching,  sorting,  extraction,  etc.  The 
new  format  enables  fabrication-dependent  material  properties  (coming  from  different 
semiconductor  fabs),  corresponding  to  the  actual  MEMS  device  being  simulated.  See 
detailed  description  below. 

-  GT  participated  in  collecting  and  providing  the  input  of  fabrication-dependent  data  into  the 

database.  The  material  data  are  being  placed  in  the  Material  Property  Manager  developed  by 
CFDRC  in  this  task. 

Property  Manager  (PM)  is  a  database  system  which  will  help  the  user  to  manage  various 
materials  and  their  associated  properties.  The  basic  database  functions  include  Add,  Modify, 
Delete,  and  Sort.  The  stand  alone  version  of  PM  provides  the  user  a  simple  and  fast  way  to 
access  the  material  database.  Also,  it  can  be  integrated  to  CFD-GUI  and  other  software  packages 
in  order  to  define  materials  which  are  part  of  a  problem  setup.  All  data  in  the  Property  Manager 
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will  be  treated  as  generic  data.  No  physical  meaning  is  programmed  into  PM  so  the  database  can 
be  expanded  virtually  unlimited. 


Data  Structure 

The  data  structure  of  the  material  database  can  be  represented  as  the  following  figure. 


PM 


Material  1 


Attributes 

|->Key  1 - ►Value 

L>Keyn  - ►Value 

Properties 


Property  1 

Evaluation  Method  1 
Variable  1 


I — ►Variable  n 


L 


Evaluation  Method  n 


u 


Property  n 


Material  n 


The  graphical  user  interface  of  PM  reflects  the  above  structure. 

Usage 

Using  PM  is  similar  to  using  a  word  processor.  The  user  can  use  the  File  menu  to  open,  save  and 
import  database  files.  The  title  of  the  window  shows  the  current  working  file  name  and  the  status 
of  that  file.  There  is  New  function  to  clear  all  materials  and  there  is  a  warning  for  unsaved 
database. 

Under  Edit  menu,  the  user  can  access  Add  and  Delete  functions.  When  the  user  selects  an  item, 
he  can  edit  the  name  and  value  in  the  right  side  edit  panel. 

The  Preference  function  provides  two  sorting  methods.  The  default  is  Sort  by  Name.  The  user 
can  also  sort  the  materials  by  the  Attribute  Type.  This  option  is  only  available  when  there  is  at 
least  one  Type  in  the  database. 


The  Material  node  and  Property  node  only  contain  its  name.  The  Attribute  node  was  constituted 
by  Key  and  value.  Besides  the  name,  the  user  can  specify  if  it  is  the  preferred  method  in  the 


Evaluation  Method  node.  In  the  Variable  level,  the  user  can  enter  the  variable  name,  variable 
type  (Integer,  Float  or  String)  and  its  value. 

All  of  the  above  functions  can  also  be  accessed  through  the  toolbar.  The  user  can  even  access  the 
Edit  function  by  a  right  mouse  click  to  pop  up  a  menu. 

The  Refresh  function  was  designed  for  using  after  the  user  has  added  new  items  and  wants  to 
sort  the  materials  using  the  current  sorting  method. 

All  of  the  tip  and  error  messages  are  displayed  in  the  status  bar,  and  the  Copy/Paste  function  is 
available  for  user  convenience. 

2.5  Slip  Wall  for  Rarefied  Gases 

In  the  usual  macroscopic  analysis  of  transfer  processes,  fluid  media  are  treated  as  continua,  and 
macroscopic  properties  such  as  density,  velocity  and  temperature  are  assumed  to  vary 
continuously  in  time  and  space.  However,  for  very  low  pressure  and/or  for  very  small  channel 
dimensions,  the  situation  may  be  different. 

As  an  example,  let  us  consider  air  flowing  through  a  channel  of  height  H.  If  H  is  large  compared 
to  the  mean  free  path  of  air  molecules  (at  sea  level  condition,  the  mean  free  path  of  air,  X,  is 

about  6T0~^  m),  the  continuum  theory  postulates  existence  of  boundary  layer  at  the  channel 
wall,  i.e.,  the  velocity  of  air  at  wall  is  the  same  as  that  of  wall.  However,  if  the  channel  height  is 
of  same  order  as  the  intermolecular  collisions  become  less  frequent  and  molecules  arriving  at 
the  wall  may  not  reach  equilibrium  with  the  wall.  The  velocity  of  air  at  wall  may  be 
discontinuous  and  the  “no  slip”  wall  condition  is  no  longer  valid.  The  so-called  “slip”  boundary 
condition  must  be  applied. 

The  non-dimensional  parameter,  K,„  Knudsen  number,  can  be  used  to  describe  such  situation. 
The  definition  of  the  Knudsen  number  is 


K 


n 


X 

L 


(2.4) 


where,  L  is  the  system  characteristic  length.  Based  on  the  Knudsen  number,  fluid  flow  can  be 
distinguished  into  different  regimes: 


Continuum  regime: 
Slip  flow  regime: 
Transitional  regime; 
Free,  molecular  flow: 


Kn  <  0.01 

0.01<K„<0.1 
0.1  <K„<3 
K„>3 


This  description  focuses  on  the  slip  flow  regime.  Following  [Rohsenow  and  Choi,  1961],  we 
assume  that  the  gas  velocity  at  wall  is  Vg  ,  and  the  wall  velocity  is  Vw  The  shear  stress  of  flow  at 
the  wall  is 
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=  (2-5) 

where,  ^  is  molecular  viscosity  of  gas,  3  is  proportionality  constant.  On  the  other  hand,  from  the 
theory  of  molecular  dynamics,  molecules  near  the  wall  have  the  last  collision  with  wall  at  an 
2A 

average  distance  of  — .  The  number  of  molecules  striking  the  wall  per  unit  area  and  time  is 
Vn  — 

— .  Where,  n  is  number  density  of  molecules  and  V  is  average  velocity  of  a  molecule.  Thus, 
4 

the  tangential  momentum  of  molecules  when  they  hit  wall  is 


m  Vn  + 


2X  8u 
3  3y 


where,  m  is  the  mass  of  gas  molecule.  If  the  molecules  striking  the  wall  transfer  a  fraction  of  F 
tangential  momentum,  we  can  balance  the  difference  of  momentum  with  wall  shear  stress: 


_mnV(_-  2Xdu  ]  du 
F— —  Vo+— — -V^ 

41  3  dy  j  dy 


Again,  from  molecular  dynamics  theory. 


p  =  nmV- 
3 


Therefore,  we  have 


A  more  refined  analysis  gives 


V  V  _22-F,au 


("2-F^2 


3  V  F  ;3 


(2.10) 


3“  F 


(2.11) 


The  values  of  F  for  different  gas-surface  pairs  are  listed  in  the  following  table. 
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Table  2.1  Values  of  F  for  different  gas-surface  pairs. 


Gas 

Surface 

F  (Accommodation  Coefficient) 

Air 

Oil 

0.90 

He 

Oil 

0.87 

H2 

Oil 

0.93 

CO2 

Oil 

0.92 

Air 

Hg 

1.00 

Air 

Machined  Brass 

1.00 

He 

Polished  Ag20 

1.00 

H2 

Polished  Ag20 

1.00 

O2 

Polished  Ag20 

0.99 

Air 

Polished  Ag20 

0.98 

The  gas  slip  velocity  at  wall  is  given  by 


(2.12) 


The  equation  above  can  be  used  with  standard  molecular  viscosity  value  in  the  slip-flow  regime 
as  mentioned  above,  that  is  for  Kn  <  0.1.  For  some  particular  cases  of  higher  Knudsen  numbers 
and  specific  geometries,  the  user  may  modify  the  viscosity  value  in  CFD-ACE+,  to  so-called 
effective  viscosity,  using  the  built-in  formula: 


lieff  = 


It 


(i  +  skS) 


(2.13) 


where  //  is  the  default  viscosity  value,  a  and  b  are  constants  given  by  the  user. 

For  example,  in  the  paper  [Veijola,  1995],  an  effective  viscosity  model  for  the  Poiseuille  flow 
(squeeze-film  damping)  has  been  derived,  with  the  constants  a  -  9.638,  Z?  =  1.159.  The  model  is 
valid  for  Knudsen  numbers  in  the  range  0  <  Kn  <  880,  with  less  than  5%  difference  from  the 
solution  of  Boltzmann  transport  equation  for  that  case.  For  a  Couette-type  flow,  that  is,  a  plate 
moving  horizontally  over  a  flat  surface,  with  a  very  small  gap  between  the  plate  and  the  surface, 
the  other  set  of  constants  has  been  proposed  in  [Veijola,  2000]:  a  =  2,  6  =  1. 

This  slip  wall  model  has  been  incorporated  into  CFD-ACE+  and  tested  on  several  cases.  As  an 
example,  the  flow  of  Helium  in  a  microchannel  was  analyzed.  The  channel  had  a  rectangular 
cross-sections  with  dimensions  7500  *  52.25  *  1.33  pm.  The  exit  pressure  was  held  constant  at 
108000  Pa,  while  the  upstream  or  inlet  pressure  was  varied.  The  experimental  data  for  mass  flow 
rate  for  He  were  published  in  [Arklic,  1994]  for  inlet  to  outlet  pressure  ratio  (Pinie/Pexu)  of  1.1  to 
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2.6.  For  the  flow  conditions  and  the  microchannel  dimensions,  the  Knudsen  number  is  about 
0.08,  i.e.,  in  the  slip  region.  A  3-D  grid  with  60*10*5  cells  was  built  for  the  channel.  Two  sets  of 
results  were  generated.  One  for  “no  slip”  wall  condition  and  another  for  “slip”  wall  condition 
using  the  model  described  above.  The  computed  mass  flows  with  the  both  wall  conditions  are 
plotted  in  Figure  9,  together  with  experimental  data  from  MIT  [Arklic,  1994].  As  can  be  seen, 
the  no-slip  wall  condition  generates  higher  wall  shear  with  a  consequent  lower  mass  flow  rate  at 
a  given  pressure  ratio.  With  the  slip  wall  condition,  the  shear  at  wall  is  lower  and  the  mass  flow 
rate  increases.  It  can  be  seen  in  Figure  9  that  the  numerical  results  with  slip  wall  condition  are 
very  close  to  the  experimental  data  for  this  case. 


Figure  9.  Comparison  of  mass  flow  rfate  for  Kn  =  0.08. 

2.6  Validation  of  CFD-ACE+  for  High-Fidelitv  Simulations 

The  single  synthetic  jet  measurements  performed  at  GT  were  used  to  compare  with  the  2D/3-D 
simulations  of  these  devices,  performed  at  CFDRC  using  CFD-ACE+.  First,  a  single  synthetic  jet 
was  used  which  was  characterized  in  great  detail  at  GT.  The  geometry  and  boundary  conditions 
used  for  simulations  were  provided  by  the  GT  team.  The  diameter  of  the  orifice  was  6.35  mm 
and  the  diameter  of  the  actuator  at  the  bottom  19  mm.  The  height  of  the  actuator  was  95  mm. 
The  oscillation  frequency  of  the  bottom  wall,  diaphragm,  was  930  Hz,  700  Hz,  and  450  Hz. 

The  dependence  of  the  actuator  stroke  length  LJD  on  the  (rms)  cavity  pressure  for  actuation 

frequency  930  Hz  is  shown  in  Figure  10.  One  can  see  that  the  agreement  between  the  measured 
data  and  simulation  results  is  very  good.  The  only  difference  is  a  kind  of  discontinuity  in  the 
experimental  results  at  the  pressure  value  about  0.1  psig.  This  effect  is  probably  caused  by 
experimental  setup  details  rather  than  by  basic  physical  behavior  of  the  synthetic  jet. 
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Figure  10.  Actuator  stroke  length  L„/D  vs.  the  (rms)  cavity  pressure,  at  frequency  930  Hz, 

measured  with  PIV  technique  at  GeorgiaTech,  and  obtained  from  CFD-ACE+  3-D  simulations  at 

CFDRC. 

Measured  and  simulated  results  of  vortex  formation  are  shown  in  the  form  of  contour  maps  in 
Figure  11.  The  measured  ensemble-averaged  vorticity  (w)  at  t*  =0.5,  for  L^jD  values  between 
0.4  and  2.2,  at  the  two  formation  frequencies  are  presented,  while  the  simulation  results  are  for 
the  frequency  fi  =  930  Hz  and  four  selected  L^jD  values.  Both  vertical  and  horizontal  axes  on 
both  sets  of  results  have  the  same  range.  This  time,  with  the  correct  value  of  the  orifice  thickness, 
the  agreement  between  PIV  and  CFD  results  is  very  good. 
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(a)  PIV  Measurements  (b)  CFD  Simulations 


Figure  11.  Contour  maps  of  phase  averaged  dimensional  azimuthal  vorticity  at  t*  =  0.5. 
(a)  PIV  measurements  (fi  =  930  Hz,  f2  =  450  Hz),  (b)  CFD-ACE+  simulations. 

3  INTERFACES  TO  CAD/EDA  TOOLS 

3.1  CFD-Micromesh:  Automatic  Generation  of  3-D  Device  Models  from  Layouts 


This  project  demonstrated  the  feasibility  of  automatic  generating  3-D  MEMS  structure  models 
directly  from  layouts  imported  from  commercial  EDA  tools.  A  new  software,  called  CFD- 
Micromesh,  was  developed  during  this  project  at  CFDRC,  for  automatic  3-D  solid  model 
generation  and  meshing  from  MEMS  layouts  imported  directly  from  EDA  systems.  The  popular 
layout  description  formats,  CIF  and  GDSII,  are  imported  easily  into  CFD-Micromesh,  allowing 
coupling  the  software  with  several  commercial  EDA  design  tools.  The  new  tool  is  fast  and 
equipped  with  user-friendly  graphical  user  interface  (GUI)  -  see  Figure  12. 
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In  Figure  12,  we  present  a  layout  of  a  MEMS  comb  drive  resonator,  imported  from  CIF  or 
GDSII  file.  To  assign  physical  meaning  to  CIF  or  GDSII  layers,  a  technology  file  is  used  in  a 
special  CFD-Micromesh  format,  for  example  MUMPs.umt  file  for  the  MUMPs  technology. 


Figure  12.  Layout  of  a  comb  resonator,  imported  from  CIF  file  into  CFD-Micromesh. 


By  clicking  only  a  single  button  in  CFD-Micromesh,  a  three-dimensional  solid  model  is 
generated  fully  automatically  from  layout.  Building  the  3-D  model  of  the  comb-drive  resonator, 
presented  in  Figure  13,  took  only  40  seconds  on  a  700  MHz  Pentium  PC.  A  “3-D  Quick  View” 
option  in  CFD-Micromesh  builds  3-D  models  of  planar  structures  in  1-2  seconds. 
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Figure  13.  A  3-D  solid  model  generated  automatically  from  layout  in  Figure  12. 

Sixty  more  seconds  were  needed  to  generate  automatically  a  full  3-D  unstructured  computational 
mesh  in  DTF  format,  again  only  by  clicking  a  single  button  in  CFD-Micromesh.  Such  a  mesh  for 
the  analyzed  comb-drive  resonator  is  shown  in  Figure  14. 


Figure  14.  A  3-D  computational  mesh,  generated  automatically  from  layout  in  Figure  12. 
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The  automatically  generated  three-dimensional  device  model  with  the  unstructured  3-D 
computational  mesh  in  DTF  format  can  be  imported  directly  into  CFDRC's  simulator  CFD- 
ACE+  for  subsequent  high-fidelity  3-D  simulations,  both  steady-state  or  transient. 

CFD-Micromesh  was  recently  extended  with  a  capability  to  read  a  fabrication  process  (layer 
structures)  description  in  a  standard  format  called  SIPPs  (Standard  Interconnect  Performance 
Parameters).  The  format  is  being  developed  by  several  big  semiconductor  companies,  like  Si2, 
AMD,  Motorola,  IBM,  Cadence,  Avanti,  Mentor  Graphics,  and  others,  with  participation  of 
CFDRC  (more  info:  httD://www.si2.org/sipps/  ).  SIPPs  layer  data  are  typically  imported  in 
connection  with  a  layout  in  GDSII  format,  which  makes  together  a  complete  3-D  description  of  a 
design.  This  capability  makes  CFD-Micromesh  a  very  attractive  modeling  tool  also  for  VLSI 
electronics  industry,  in  particular  for  analysis  of  interconnects  and  PCBs. 

The  new  version  of  CFD-Micromesh  imports  CIF  and  GDSII  layouts,  reads  the  SIPPs  process 
description  format,  and  includes  direct  interface  to  Virtuoso^^  layout  editor  from  Cadence.  The 
evaluation  version  (available  from  http://www.cfdrc.com/~micromesh)  builds  automatically  3-D 
solid  models  from  layout  and  process  data,  and  allows  the  user  to  view  them  interactively. 

Finally,  CFDRC  has  negotiated  an  agreement  with  Cadence  Design  Systems,  Inc.  to  develop 
interfaces  and  integration  of  CFDRC's  MEMS  simulation  software  with  Cadence  CAD  tools. 
CFD-Micromesh  will  be  a  part  of  Cadence  CAD  Environment,  used  as  a  “3-D  Project  Viewer” 
from  Virtuoso  layout  editor  from  Cadence  (see  Figure  5). 

The  “3-D  Project  Viewer”  interface  to  Cadence  has  been  included  into  the  commercially 
available  CFD-Micromesh  package.  The  appropriate  linking  procedures  (in  Cadence’s  SKILL 
Script)  and  basic  documentation  of  the  interface  are  available  in  the  CFD-Micromesh  package. 

The  interface  and  all  of  the  involved  procedures  of  coupling  CFD-ACE-i-  with  SABER  system- 
level  simulator  have  been  described  in  detail  in  the  following  publications: 

•  M.  Furmanczyk,  M.  Turowski,  E.  Yu,  and  A.  Przekwas,  "Coupled  System-Level  and 
Physical-Level  Simulation",  ASSURE  2000  (Avant!  Saber  Simulator  Users  Resource 
Meeting),  Portland,  Oregon,  May  10-12,  2000. 

•  M.  Furmanczyk,  M.  Turowski,  E.  Yu,  and  A.  Przekwas,  "Coupled  System-Level  and 
Physical-Level  Simulation",  Int.  Conf.  MSM  2001,  Hilton  Head  Island,  South  Carolina, 
USA,  19-21  March  2001. 

3.2  Schematic  Driven  Design  Entry  for  MEMS 

The  goal  of  this  task,  performed  mainly  by  the  CMU  team  in  collaboration  with  CFDRC,  was  to 
develop  an  integrated  CAD  environment  composed  of  schematic  driven  design  entry  and 
automated  generation  of  3-D  mesh  and  solids  description  for  rapid  design  of  MEMS  structures. 

Design  process  of  MEMS  structures  is  typically  very  time-consuming,  especially  in  terms  of 
design  verification  and  redesign  issues.  To  make  it  faster  and  more  effective,  we  have  proposed 
to  introduce  a  schematic  entry  into  the  existing  design  flow  sequence,  similarly  to  the  established 
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EDA  procedures  for  VI, SI  integrated  circuits.  The  modified  MEMS  design  flow  would  consist 
of  tiic  following  steps: 

1)  introducing  the  idea  in  a  form  of  MEMS  schematic, 

2)  verification  by  means  of  complex  3-D  and  multidomain  simulations,  (optionally  go  to 
step  1 ) 

3)  layout  generation  and  verification,  (optionally  go  to  step  I). 

In  the  frame  of  this  project,  we  have  developed  a  CAD  environment  allowing  to  generate  3-D 
mesh  and  solids  description  in  PLOT3-D  and  DTE  formats  (viewed  in  CFD-VIEW  package) 
from  MEMS  schematics  (created  in  Cadence  Schematic  Editor). 


CAF^ENC'i:; 

(‘DvinujnTftMir 


procedures 


Bxiemal  .Ml  fl  ware 
C'i-O-.ACE 


3-D  simulation 


Figure  15.  The  full  design  flow  for  MEMS  structures.  The  shaded  areas  indicate  the  field  of  interest 

in  this  project. 
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CMU  NODAS  CAD  tools  allow  fast  design  and  verification  of  MEMS  devices  built  from  a  set 
of  “atomic”  elements  like  beams,  plates,  combs  and  anchors.  Figure  15  describes  in  general  view 
the  design  flow  with  the  shaded  areas  indicating  the  field  of  interest  in  the  frame  of  this  part  of 
the  project.  The  design  of  MEMS  structures,  where  the  influence  of  the  surrounding  environment 
on  device  functioning  has  to  be  taken  into  account,  requires  3-D  and  multidomain  simulations.  In 
the  case  of  e.g.  CMOS  ultrasonic  structures,  air  damping  effect  strongly  determines  the  device 
performance. 

The  developed  in  Cadence  SKILL  language  CAD  procedures  use  the  description  of  NODAS 
symbols  parameters  taken  from  the  schematic  database  separately  for  each  type  of  symbols. 
Then,  a  3-D  rectangular  mesh  description  is  created  separately  for  each  set  of  symbols  as  well. 
The  output  format  can  be  PLOT3-D  (structured  grid  with  blanking)  or  DTF  format  developed  at 
CFDRC.  Any  changes  of  3-D  description  are  made  by  editing  parameter  values  in  Edit 
Properties  window  invoked  from  Cadence  Schematic  Editor  window.  Redesign  of  such  a 
structure  in  terms  of  changes  in  geometric  dimensions  takes  up  to  several  seconds  only.  The 
developed  CAD  procedures  are  accessed  by  the  user  through  the  GUI  in  Cadence  environment. 
An  exemplary  design  views  of  folded  flexure  resonator  in  MUMPS  technology,  as  well  as  the 
GUI  elements  are  shown  in  Figure  16. 


Figure  16.  Schematic,  layout  and  3-D  visualization  views  of  MUMPS  folded  flexure  resonator. 
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As  a  demonstration,  the  new  technology  has  been  applied  to: 

•  3-D  simulations  of  CMOS  MEMS  ultrasonic  sensor. 

•  A  folded  beam  MEMS  structure 

3.3  CAD  Tools  Integration  for  MEMS  Design 

This  sub-task  aimed  at  developing  an  integrated  CAD  environment  that  combines  the  ease  of 
schematic  driven  design  entry  and  automated  generation  of  3-D  mesh  and  solids  for  accurate 
continuum  simulation  of  MEMS  devices.  This  effort  links  the  NODAS  schematic  simulation 
environment  with  the  3-D  visualization  and  device  verification  capabilities  in  CFD-ACE-+. 

There  are  two  focus  points  for  this  sub-task: 

•  3-D  visualization  in  which  the  schematic  needs  to  be  translated  into  a  meshed  geometry 
for  visualization  in  CFD-VIEW. 

•  Device  verification  in  which  both  the  schematic  geometry  as  well  as  the  source 
excitations  on  the  schematic  need  to  be  translated.  The  source  excitations  in  the 
schematic  become  boundary  conditions  for  the  CFD-ACE+  solver. 

We  have  developed  a  new  data  format  (called  GBV  -  Geometry,  Boundary  and  Volume  conditions) 
for  interchange  of  simulation  data  between  CADENCE  Design  Framework  II  and  CFDRC  ACE 
simulator.  Based  on  this  format.  Cadence  design  data  extraction  procedures  (SKILL)  and  program 
(C)  translating  data  from  GBV  format  to  CFDRC-DTF  format  have  been  created.  The  design  flow 
is  based  on  generation  of  text  GBV  file  directly  from  CADENCE  Schematic  Editor  by  reading  out 
schematic  database,  running  layout  generation  and  meshing  procedures.  The  layout  meshing  is 
accomplished  using  the  algorithms  detailed  in  [Hasnain  Lakdawala,  Bikram  Baidya,  Tamal 
Mukherjee  and  Gary  K.  Fedder,  “Intelligent  Automatic  Meshing  Of  Multilayer  CMOS 
Micromachined  Structures  For  Finite  Element  Analysis,”  in  Proc.  MSM  1999,  pp.  297-300].  The 
geometry  data  obtained  after  layout  meshing  is  combined  with  boundary  and  volume  conditions 
data  taken  from  schematic  view  in  order  to  create  a  file  with  a  data  description  in  GBV  format. 
Next,  the  text  file  (GBV)  is  translated  into  binary  DTF  file,  which  is  required  by  the  CFD-ACE 
solver.  Generation  of  the  new  DTF  file  can  be  performed  by  changing  values  of  symbol  parameters 
in  the  Cadence  Schematic  Editor.  GUI  procedures,  integrating  CADENCE  tools  with  the  CFDRC 
ACE  simulator  for  the  mentioned  design  flow,  have  been  also  developed. 

The  examples  in  the  figures  below  demonstrate  the  capability  of  the  resulting  link  between 
NODAS  and  CFD-ACE-i-.  Figures  17  and  18  demonstrate  the  visualization  capabilities.  Figures 
19-21  demonstrate  the  capability  of  passing  schematic  sources  as  boundary  conditions  to  CFD- 
ACE+. 
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Figure  18.  CFD-GUI  view  of  the  beams  and  beam  corners. 
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Figure  19.  Snapshots  of  Cadence  Schematic  Editor  with  anchor-beam-beam- 
dc_force_source  test  design  and  its  view  in  CFDRC-GUI  after  translation  to  DTF  format. 


Figure  20.  Snapshot  of  CFDRC-GUI  with  test  design  showing  its  layers  (eight  CMOS 
layers)  and  Boundary  Condition  description  for  the  wall  connected  to  the  anchor  on 
schematic  view  (Fixed  Displacement  dx=0,  dy=0,  dz=0). 
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Figure  21.  Simulation  results  for  the  test  design:  analysis  type:  steady,  dc_force  value: 

lOOuN  new,  reduced  and  compact  models 

4  NEW  REDUCED  AND  COMPACT  MODELS 

4.1  Air  Damning  Models 

4.1.1  Nonlinear  Compact  Model  of  Squeeze  Film  Damping 

MEMS  devices  are  characterized  by  very  small  gaps  between  the  moving  elements  and  the  fixed 
parts.  A  gas  (air)  film  between  two  closely  spaced  plates  is  squeezed  and  produces  forces  that 
oppose  the  motion  of  the  plate.  The  viscous  and  compressibility  effects  dissipate  the  energy  of 
the  moving  plate,  which  is  known  as  squeeze-film  damping  (Figure  22a).  Several  authors  have 
already  analyzed  the  squeeze-film  behavior,  but  all  the  solutions  have  been  shown  only  for  small 
amplitudes  of  the  motion.  An  analytical  solution  of  the  linearized  compressible  isothermal 
Reynolds  equation  was  proposed  by  Blech  in  [Blech,  1983],  under  the  assumption  that  the 
motion  of  the  plates  and  the  pressure  variations  are  small.  The  analytical  solution  for  the 
squeeze-film  force  consists  of  two  components:  one  is  in  phase  with  the  plate  velocity  (the 
viscous  damping  force),  the  other  is  in  phase  with  the  plate  displacement  (the  spring  force  due  to 
the  compressibility  of  air). 
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(a)  (b) 

Figure  22.  Squeeze  film  between  plates  (a),  and  its  equivalent-circuit  model  (b)  The  arrows  show 

the  role  of  CFD-ACE+  high-lldelity  simulations. 


On  the  basis  of  the  Blech's  solution,  an  equivalent-circuit  model  of  squeeze-film  damping  was 
proposed  in  [Veijola  1995],  realized  with  a  ladder  of  RL  (resistor  and  inductor)  branches  (Figure 
22b)  satisfying  the  frequency  dependence  of  both  analytical  components:  damping  force  and 
spring  force.  This  model  was  also  derived  with  the  assumption  of  small  amplitudes  of  motion 
and  small  pressure  changes,  and  for  such  conditions  a  good  agreement  of  the  circuit  model  with 
measurements  was  demonstrated  in  [Veijola  1995].  A  modification  of  the  R  and  L  values  in  the 
equivalent  circuit,  as  nonlinear  functions  of  displacement,  was  also  proposed  in  [Veijola  1995], 
but  neither  was  it  verified  nor  any  large-amplitude  results  were  shown. 

The  same  problem  was  investigated  by  Yang  and  Senturia  in  [Yang  1996],  where  squeeze-film 
effects  were  modeled  numerically  under  small-amplitude  sinusoidal  conditions  by  a  finite 
element  package,  as  well  as  under  large  amplitude  motion  by  a  finite  difference  codes.  For  the 
large  displacement,  the  maximum  oscillation  amplitude  in  that  work  was  equal  to  30%  of  the 
nominal  gap  thickness,  for  which  only  numerical  finite-difference  solution  was  shown.  Also  the 
most  recent  papers  reporting  compact  analytical  models  of  squeeze  film  damping,  like  [Yang 
1997],  [Darling  1997]  ,  [Burgeois  1997]or  [Pan  1998],  are  limited  to  small  amplitude  motion  as 
one  of  their  main  assumptions,  which  may  be  not  necessarily  true  in  real  MEMS  behavior. 

In  this  project,  we  performed,  demonstrated,  and  published  simulations  of  isothermal 
compressible  squeeze-film  damping  effects  between  plates  of  MEMS  devices  for  arbitrary  large 
amplitudes,  using  a  behavioral  model  based  on  equivalent  nonlinear  RL  circuit.  The  results  of  the 
compact  model  were  verified  using  a  3-D  high-fidelity  transient  solution  of  full  Navier-Stokes 
equations,  obtained  with  CFD-ACE+  from  CFDRC.  The  compact  model,  based  on  the  format 
proposed  in  [Veijola  1995],  is  shown  in  Figure  22b.  The  component  values  in  the  RL  ladder  are: 
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where:  I  and  w  are  the  length  and  width  of  the  plate,  respectively,  g  is  the  nominal  gap  height 
(Figure22a),  z  is  plate  displacement,  and  Pa  is  the  static  pressure.  The  effective  viscosity 

in  Equation  (4.2)  depends  on  the  gap  height  g  and  the  mean  free  path  A  of  the  gas,  through  the 
Knudsen  number  =  A/g.  An  expression  for  the  effective  viscosity  is  given  in  [Veijola  1995]. 

The  coefficients  m,  n  are  odd  indices  1,3,5,...,  reflecting  the  series  expansion  form  of  the 
squeeze  film  solution.  In  practice,  only  a  few  sections  are  needed  for  sufficient  accuracy.  In  our 
simulations  we  have  used  the  first  six  RL  sections,  for  {m,n)  =  (1,1),  (1,3).  (1,5),  (3,1),  (3,3), 
(5,1). 

The  above  mentioned  electrical  equivalences  for  mechanical  and  fluidic  domains  are  used  here: 
the  voltage/current  pair  represents  velocity/force.  Hence,  the  plate  velocity  V  and  the  squeeze 
force  Fsq  in  Figure  22a  are  represented  by  the  voltage  and  the  current  ig  in  Figure  22b, 

respectively.  The  plate  displacement  z  can  be  calculated  from  the  flux  ifr  of  the  resonator 
inductance  L^pr  ,  if  we  add  to  the  equivalent  circuit  an  additional  inductor  L^pj-  =  M, 
representing  the  device  elastic  spring  constant  k.  If  the  inductor  L^pf  was  connected  between 
pins  A  and  B  in  Figure  22b,  the  displaeement  would  be:  z  =  (^  =  Lspyiij  where  denotes  the 
current  through  the  inductor  L^pj- .  If  there  is  no  spring  constant  representation  in  the  circuit,  the 
displacement  z  can  be  calculated  by  integrating  the  velocity  which  is  directly  represented  by  the 
voltage  in  the  electric  equivalent  circuit. 

SPICE  implementation  and  comparison  with  3-D  numerical  simulations 

The  equivalent  RL  model  was  implemented  in  our  case  first  in  PSpice  running  on  PC.  The 
SPICE  versions  known  to  us  do  not  have  nonlinear  models  of  resistors  R  and  inductors  L  that 
would  allow  for  direct  implementation  of  Equations  (4.1)  and  (4.2).  As  explained  above,  the  z 
value  in  the  equations  for  L  and  R  comes  either  from  current  in  another  branch,  or  from 
integrating  the  voltage  in  the  circuit.  Therefore,  in  the  SPICE  model,  the  nonlinear  elements 
and  have  to  be  represented  by  nonlinear  controlled  sources.  In  our  implementation  for 
PSpice,  we  have  used  voltage-controlled  voltage  sources  (VCVS)  for  inductors  L^„,  and 
voltage-controlled  current  sources  (VCCS)  for  resistors  „ 

To  verify  the  results  of  the  equivalent  circuit  model  for  large  amplitudes,  comparisons  between 
the  3-D  high-fidelity  simulations  and  compact  model  results  were  done,  in  frequency  and  time 
domains,  for  periodic  plate  displacements.  The  3-D  numerical  calculations  were  done  using 
CFD-ACE-I-.  They  involved  solution  of  full  Navier-Stokes  equations,  as  described  in  the 
previous  section.  Our  effort  in  this  project  demonstrated  for  the  first  time  the  squeeze  film 
behavior  in  MEMS  for  very  large  amplitudes,  up  to  90%  of  the  nominal  gap  thickness,  and 
our  results  have  been  published  in  several  papers  (see  “Publications  Resulting  from  This 
Project”). 

Example  results  obtained  with  full  3-D  Navier-Stokes  solution  and  equivalent-circuit  RL  model 
(solved  by  SPICE),  for  periodic  oscillations  of  frequency  100  Hz,  are  presented  in  Figure  23  and 
Figure  24.  The  nonlinear  compact-model  results  agree  very  well  with  3-D  results  even  for  the 
very  large  amplitudes  of  plate  motion,  accompanied  by  significant  changes  of  pressure.  In  the 
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figures,  the  displacement  amplitude  equals  90%  of  the  nominal  gap  height  (4  micrometers  in  this 
case),  and  the  resulting  relative  pressure  change  with  respect  to  static  ambient  pressure,  due  to 
the  squeeze  force,  reaches  the  maximum  of  about  3,  which  is  not  even  close  to  the  small 
amplitude  and  small  pressure  change  assumptions  made  previously. 


CFD-ACE:  3-D  model 


SPICE:  RL  model 
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Figure  23.  Relative  pressure  change  with  respect  to  static  ambient  pressure  as  a  function  of  plate 
displacement,  obtained  with  3-D  Navier-Stokes  solution  (using  CFD-ACE)  and  equivalent-circuit 
model  (solved  by  SPICE),  for  the  displacement  amplitude  equal  to  90%  of  the  nominal  gap  height. 
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Figure  24.  Relative  pressure  change  with  respect  to  static  ambient  pressure  as  a  function  of  plate 
velocity,  obtained  with  3-D  numerical  solution  and  equivalent-circuit  model,  for  the  displacement 

amplitude  equal  to  90%  of  the  nominal  gap  height. 

In  Figure  25,  the  calculated  relative  pressure  change  due  to  squeeze  effect  is  presented  for  very 
high  operating  frequency,  equal  100  kHz.  The  displacement  amplitude  is  again  90%  of  the 
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nominal  gap  height  (4  micrometers  in  this  case  as  well).  The  resulting  relative  pressure  change  is 
10  to  15  times  bigger  than  the  static  ambient  pressure,  which  means  even  bigger  changes  of 
pressure  than  for  low  fnequencies.  The  agreement  between  RL-circuit  and  3-D  model  for  the 
high  frequency  is  not  as  good  as  for  low  frequencies,  but  the  shape  of  characteristics  is  preserved 
and  we  believe  the  equiv  alent-circuit  RL  model  still  can  be  used  even  for  high  frequencies,  after 
some  parameter  fitting. 
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Figure  25.  Relative  pressure  change  with  respect  to  static  pressure  as  a  function  of  plate 
displacement  for  frequency  100  kHz  and  displacement  amplitude  90%  of  the  nominal  gap. 


S»aber-MAST  implementation 


The  equivalent-circuit  model  of  the  squeeze  film  behavior  has  been  also  implemented  for  the 
Saber  simulator  with  models  written  in  MAST,  an  analog  hardware  description  language 
[MAST,  1997].  The  Saber-MAST  version  of  the  behavioral  squeeze  film  model  uses  directly 
mechanical  quantities:  position  as  input  (across  variable)  and  damping  force  as  output  (through 
variable),  which  is  compatible  with  hierarchical  representation  for  system-level  design  of 
MEMS,  like  the  NODAS  methodology  developed  at  Carnegie  Mellon  University  [Vandemeer, 
1998],  partly  also  in  the  frame  of  the  DARPA  Composite  CAD  Program  .  The  developed  model, 
together  with  the  simple  test  circuit  containing  also  a  mechanical  translation  source  available  in 
Saber,  is  shown  in  Figure  26.  The  description  of  the  model  is  included  in  a  MAST  template 
called  squeezelD.  This  version  of  model  is  only  one-dimensional  (ID),  or  for  one  degree-of- 
freedom  (DOF).  The  squeeze  damping  model  in  one  direction  has  only  two  mechanical  nodes 
(say,  "top"  and  "bottom"),  with  their  displacements  (x_t,  x_b)  as  across  variables,  and  resulting 
squeeze  force  as  through  variable. 
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Figure  26.  Squeeze  film  behavioral  model  in  Saber  simulator,  together  with  mechanical 

displacement  source  for  testing. 

Similarly  as  in  SPICE,  the  models  of  resistors  R  and  inductors  L  available  in  Saber  libraries  do 
not  allow  for  direct  implementation  of  Equations  (4.1)  and  (4.2).  Therefore,  special  MAST 
templates  had  to  be  written  for  the  nonlinear  elements  „  and  Rm,n  ^he  equivalent  circuit. 

The  z  value  denoting  plate  displacement  has  been  directly  inserted  into  the  MAST  templates, 
called  Lsq  and  Rsq,  thanks  to  the  possibility  of  mixing  mechanical  and  electrical  signals  in  one 
model.  The  MAST  code  of  the  squeezelD  template,  modeling  the  squeeze  film  behavior  with 
mechanical  inputs  (position)  and  output  (force),  is  shown  in  Figure  27.  The  template  is  generated 
automatically  by  a  special  routine  written  in  C,  taking  as  input  all  the  basic  parameters  of  the 
plate,  gap,  gas  pressure,  etc.  The  same  C  program,  after  switching  one  parameter,  generates 
automatically  input  files  for  Spice  simulations  as  well. 

Behavioral  characteristics  obtained  with  Saber  using  the  MAST  implementation  of  the  RL- 
circuit  are  of  course  almost  identical  as  the  ones  obtained  with  Spice,  presented  on  the  right-hand 
side  of  Figures  23-25.  As  another  example,  in  Figure  28,  we  show  transient  waveforms  of  the 
squeeze-film  force  in  response  to  oscillatory  movement  of  the  top  plate  with  amplitude  equal 
60%  of  the  gap  height  (amplitude  =  2.4  |xm),  calculated  both  with  3-D  CFD-ACE  model  and  the 
Saber  behavioral  model. 
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###  MAST  code  for  RL  Squeeze -Film  Model:  squeezelD . sin 
# 

#  Marek  Turowski,  (C)  1998,  CFD  Research  Corporation,  Huntsville,  AL 

# 

#  MOVING  PLATE: 

#  gap  =  4.00e-006  m,  length  =  2.96e-003  m,  width  =  1.78e-003  m 

#  Pgas  =  430.0  Pa,  Area  =  5.27e-006  m2,  Static_Force  (Pgas*A)  =  2.27e-003  N 

element  template  squeezelD  x__t,  x_b  =  gap 

translational_j)os  x_t,  x_b  #  displacement  of  ’top’  and  'bottom'  plate 

number  gap  =  4u  #  gap  between  plates  [m]  (with  default  value) 

{  #  start  template  body 

branch  Fsqz_x  =  frc_N (x_t->x_b)  #  through  variable:  ID  squeeze  force  [N] 

val  nu  F_rel  #  Relative  Squeeze-Force,  F/Fa  =  Fsq/ (Area*Pamb) 

F_rel  =  i (dxdt.v) /2 .265584e-003 

#  dxdt  -  template  converting  mechanical  movement  to  voltage 

dxdt.v  x_t  x_b  vel  0  #  vel  =  velocity  of  plates  =  voltage  for  RL  model 

Lsq.Ll_l  vel  4  x_t  x_b  =  Llin  =  6 , 717990e+002 

Rsq.Rl_l  4  0  x_t  x_b  =  Rlin  =  2 . 286073e+017 

Lsq.Ll_3  vel  5  x_t  x_b  =  Llin  =  6 . 046191e+003 

Rsq.Rl_3  5  0  x_t  x_b  =  Rlin  =  6 . 428883e+018 

Lsq.Ll__5  vel  6  x_t  x_b  =  Llin  =  1 . 679497e+004 

Rsq.Rl_5  6  0  x_t  x_b  =  Rlin  =  4 . 2 14366e+019 

Lsq.L3_l  vel  7  x_t  x_b  =  Llin  =  6 . 046191e+003 

Rsq,R3__l  7  0  x_t  x_b  =  Rlin  =  1 . 41457 8e+019 

Lsq.L3_3  vel  8  x_t  x_b  =  Llin  =  5 . 441572e+004 

Rsq.R3_3  8  0  x_t  x_b  =  Rlin  =  1 . 666548e+020 

Lsq.L5_l  vel  9  x_t  x_b  =  Llin  =  1 . 679497e+004 

Rsq.R5_l  9  0  x_t  x_b  =  Rlin  =  1 . 064511e+020 

Fsqz_x  =  i (dxdt.v)  #  template  equation  for  through  variable 

}  #  end  template  body 

Figure  27.  The  MAST  code  of  the  squeezelD  template. 
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CFD-ACE:  3-D  model  Saber-MAST  model 


Figure  28.  Calculated  transient  waveforms  of  the  squeeze-film  force  in  response  to  oscillatory 
movement  of  the  top  plate  with  amplitude  equal  60%  of  the  gap  height  (amplitude  =  2.4  pm). 


The  equivalent-circuit  format  of  a  compact  (or  behavioral)  model  has  obvious  advantages:  it  can 
be  readily  implemented  in  various  system-level  simulators,  among  which  SPICE  and  Saber  seem 
to  be  the  most  popular  examples.  With  such  circuit  simulators,  linear  and  nonlinear  device 
analyses,  both  in  frequency  and  time  domains,  can  be  easily  performed.  In  addition  to  the 
standard  analyses,  advanced  circuit  simulators  offer  noise  analysis,  sensitivity  analysis,  stability 
analysis,  Monte  Carlo  simulation,  electro-thermal  simulation,  optimization,  parameter  sweeps, 
and  ready  graphical  processing  of  the  simulated  data.  When  the  MEMS  device  itself  has  an 
electrical  interface,  an  analysis  with  the  surrounding  interface  electronics,  including  feedback,  is 
easily  achieved.  Moreover,  the  analysis  is  not  restricted  to  a  single  component,  but  complex 
systems  including  MEMS  can  be  build  by  means  of  a  large  network. 

A  successful  circuit-level  modeling  of  entire  micromechanical  accelerometer,  including  the 
circuit  submodels  for  squeezed  gas  film,  was  presented  in  [Veijola  1998(a)].  The  compact 
squeeze-film  model  presented  here,  in  the  form  of  nonlinear  RL  circuit,  is  very  promising  not 
only  for  parallel  motion  of  a  rectangular  plate,  but  also  for  comb  drives  and  micromirrors.  On  the 
basis  of  the  solution  given  in  [Darling  1997],  the  RL  equivalent-circuit  model  can  be  adapted  to 
reflect  more  complex  boundary  (venting)  conditions  as  well  as  the  tilting  motion  of  a  plate  (e.g. 
for  torsional  micromirrors),  including  two  modes  of  the  motion:  normal  to  the  surface  and  a 
rotating  motion.  The  RL  compact  model  allowed  also  for  successful  modeling  of  a  plate  with 
holes,  by  fitting  parameters  for  all  the  R  and  L  elements  through  optimization.  Such  a  procedure 
was  shown  in  [Veijola  1998(b)]  using  the  AFLAC  circuit  simulator. 

The  use  of  high-fldelity  3-D  numerical  simulations,  together  with  the  procedures  described 
above,  allow  to  obtain  very  useful  compact/behavioral  models  of  squeeze-film  effects  for 
wide  range  of  MEMS  devices  and  driving  conditions. 
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4.1.2  Equivalent  Circuit  Model  of  Lateral  Viscous  Damping 


For  our  compact  model  of  the  parallel  viscous  damping,  an  analytic  expression  derived  in 
[Veijola  2000]  was  used,  which  models  the  damping  force  of  an  oscillating  plate  as  a  function  of 
frequency.  This  general  dynamic  model,  valid  in  both  the  frequency  and  the  time  domain,  is 
implemented  approximately  as  an  electrical  equivalent  circuit. 


Two  derive  the  model,  we  assumed  two  surfaces  bounding  a  flat  gas  film.  The  first  surface,  at  z 
-  d  moves  with  velocity  Vr  in  the  direction  of  the  x  axis,  while  the  second  at  z  =  0  does  not 
move,  see  Figure  29.  The  surfaces  are  assumed  large  compared  with  other  dimensions,  and  thus 
the  border  effects  are  ignored  here. 
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Figure  29.  Structure  of  the  air-gap  for  the  lateral  damping  model. 


At  low  surface  velocities  (Re«l),  the  gas  velocity  profile  can  be  assumed  to  be  linear,  and  due 
to  the  gas  rarefaction,  the  gap  width  effectively  increases  by  2X,  where  X,  is  the  mean  free  path  of 
the  gas.  According  to  [Burgdorfer,  1959],  the  shear  stress  at  one  of  the  surfaces  is 

TjA 


= 


d  +  2X 


(4.3) 


where  A  is  the  surface  area  and  //  is  the  viscosity  coefficient.  The  contribution  of  the  mean  free 
path  can  be  included  into  the  effective  viscosity: 

T] 

i^AA) 

where  Kn,  the  Knudsen  number,  is  the  measure  of  the  rarefaction  effect.  It  is  the  ratio  between 
the  mean  free  path  X  and  gap  height  d:  Kn  =  X/d.  The  resulting  damping  coefficient  ^  is  then 
simply 


^ 

V, 


(4.5) 


This  model  assumes  a  full  established  Couette  flow  ignoring  the  inertia  effect  of  the  gas.  To 
extend  the  damping  model  to  be  valid  at  higher  frequencies,  the  time-dependent  velocity  profile 
of  the  gas  must  be  considered.  The  dynamics  of  the  gas  is  modeled  with  the  one-dimensional 
diffusion  equation  [Veijola,  2000]. 

The  resulting  relation  between  the  shear  stress  x^z  and  the  velocity  Vr  in  the  frequency  domain  is 
derived  as: 
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T 


^  -^^  =  r]Aq 


cosh(^c/)  +  qAsm\\{qd) 

(1  +  q^  ^)smh{qd)  +  2qAcos\\{qd) 


(4.6) 


where  q  is  the  viscosity  coefficient,  A  is  the  surface  area,  X  is  the  mean  free  path,  q  =  -JJioJv  , 
(0  is  the  angular  frequency,  v  is  the  kinematic  viscosity  v=q/p,  and  p  is  the  gas  density.  The 
damping  coefficient  is  Re(  ^ ). 


As  proposed  in  [Veijola  2000],  the  general  dynamic  model  of  the  lateral  shear  damping,  valid  in 
both  the  frequency  and  the  time  domain,  can  be  implemented  approximately  as  an  electrical 
equivalent  circuit,  presented  in  Figure  30.  In  this  circuit,  electrical  voltage  represents  mechanical 
velocity  Vr,  and  current  represents  the  shear  stress  at  the  moving  surface,  Xxz. 
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Figure  30.  Electrical  equivalent  circuit  which  represents  the  lateral  shear  damping  forces. 


The  behavior  of  the  damping  coefficient  Re(^)  was  studied  in  [Veijola  2000]  using  the 
equivalent  circuit  model  shown  above,  as  a  function  of  frequency  when  various  parameters 
change;  such  as  the  mean  free  path  (Figure  31),  gap  height,  or  pressure.  The  equivalent-circuit 
characteristics  of  the  shear  damping  coefficient  were  compared  with  3-D  high-fidelity 
simulations  performed  at  CFDRC  using  CFD-ACE-I-.  A  series  of  3-D  and  2D  simulations  were 
run  to  extract  the  damping  coefficient  from  high-fidelity  CFD  results  (Figure  33).  To  simulate 
the  gas  rarefaction  effect  and  mean  free  path  influence  (^=^air  i*’  Figure  31),  the  slip-wall 

boundary  conditions  (BC)  were  used  in  CFD-ACE+  simulations.  Figure  32  shows  the 
corresponding  results  obtained  with  CFD-ACE+  without  slip-wall  BC  (X=0),  and  with  slip-wall 
BC  activated  (corresponding  to  ^X.aji-). 
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[  Fmax  /  Vmax  ] 


f/Hz 

Figure  31.  The  contribution  of  the  finite  mean  free  path  X  to  the  damping  coefficient  at 
atmospheric  pressure.  Gap  separation  is  1  pm  and  the  surface  area  is  1  mmxl  mm. 
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Figure  32.  Damping  coefficient  (Shear-ForceA^elocity)  characteristics  calculated  with  CFD-ACE+j 
using  Slip  Wall  boundary  condition  to  account  for  mean  free  path  contribution. 
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f  =  10  kHz,  sin.  amplitude  1  pm 


f  s  100  MHz,  sin.  amplitude  0.01  pm 


Figure  33.  Velocity  profiles  (shown  at  peak  velocity)  obtained  by  2D  simulations  with  CFD-ACE+ 
to  extract  the  damping  coefficient.  Gap  separation  is  1  pm,  the  moving  plate  edge  is  Inun. 
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The  frequency  characteristics  of  the  shear  damping  coefficient  obtained  with  2D  high-fidelity 
simulations  using  CFD-ACE+  (Figure  32)  agree  very  well  with  the  characteristics  obtained  with 
the  equivalent-circuit  (Figure  31). 

CFD-ACE+  simulations  were  also  used  to  extend  the  model  validity  for  small  plates  and  wide 
range  of  pressures.  The  gas  rarefaction  effects  were  studied  and  implemented  in  the  form  of  slip- 
wall  boundary  conditions  in  modeling.  The  equivalent-circuit  characteristics  of  the  shear 
damping  coefficient  were  compared  with  3-D  high-fidelity  simulations  performed  at  CFDRC 
using  CFD-ACE4-  for  a  wide  range  of  frequencies  and  various  ambient  pressures.  Figure  34 
shows  the  corresponding  results.  A  paper  on  this  work  was  submitted  in  June  2000  to  IEEE 
Journal  on  MEMS  (by  T.  Veijola  and  M.  Turowski). 


Figure  34.  Frequency  response  of  the  analytic  damping  coefficient  (lines),  the  real  part  of 
the  equivalent  circuit  admittance  (dashed  lines),  and  results  extracted  from  CFD-ACE+ 
simulations  ( square  symbols)  at  three  pressures.  Gap  separation  is  1  pm  and  the  surface 

area  is  1  mmxl  mm. 

The  further  work  on  enhancement  and  validation  of  the  compact/behavioral  model  of  lateral 
viscous  damping  in  MEMS  was  continued  by  the  partner  group  in  this  project  at  Carnegie 
Mellon  University  (G.  Fedder,  T.  Mukherjee,  S.  Vemuri)  -  see  below. 

4.1.3  Improvement  and  Validation  of  Viscous  Damping  Compact  Models 

Experience  with  experimental  verification  of  the  synthesis  methodologies  developed  by  the 
CMU  team  under  the  DARPA  Composite  CAD  BAA  96-16  project  on  “Foundations  for  MEMS 
Synthesis”  (MEMSYN)  demonstrated  that  the  off-the-shelf  lateral  damping  models  had  a  20% 
error  in  predicting  the  Q-factor  for  folded-flexure  resonators.  As  the  MEMSYN  project  was 
evolving  to  synthesis  of  accelerometers,  the  damping  models  governing  the  dynamics  of 
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accelerometers,  namely  squeeze  film  damping  was  considered  most  important.  We  therefore  first 
present  the  work  on  squeeze-film  modeling,  and  then  on  modeling  of  lateral  damping. 

4. 1.3.1  Parametrized  Models  of  Squeeze  Film  Damping  in  MEMS 

There  were  four  stages  to  the  development  of  parametrized  models  for  squeeze-film  damping. 
First,  we  developed  a  squeeze-film  damping  model  for  use  in  NOD  AS  as  none  had  existed  at  that 
time.  Next,  we  compared  existing  models  with  continuum  simulation  using  CFD-ACE+  to 
screen  for  model  limitations.  Next  we  developed  an  extended  model  to  overcome  the  limitations, 
and  finally,  we  verified  the  extended  model  first  using  CFD-ACE-h,  and  then  with  experimental 
measurements. 


The  starting  point  squeeze-film  damping  model  was  [J.J.Blech,  “On  Isothermal  Squeeze  Films,” 
Journal  of  Lubrication  Technology,  v.  105,  1983  ,pp.  615-620],  with  extensions  into  an 
equivalent  circuit  model  [T.  Veijola,  H.  Kuisma,  J.  Lahdenpera,  and  T.  Ryhene,  “Equivalent- 
Circuit  Model  of  Squeezed  Gas  Film  in  a  Silicon  Accelerometer,”  Sensors  and  Actuators  A,  vol. 
48,  1995,  pp.239-248]  and  to  a  non-linear  behavioral  model  [M.  Turowski,  Z.  Chen  and  A. 
Przekwas,  “Squeeze  Film  Behavior  in  MEMS  for  Large  Amplitude  Motion  -  3-D  Simulations 
and  Nonlinear  Circuit/Behavioral  Models,Proc.  lEEE/VIlJF  Inti.  Workshop  on  Behavioral 
Modeling  and  Simulation  (BMAS  ’98),  Orlando  FL,  Oct  27-28  1998]. 

As  these  models  were  derived  using  electrical  equivalent  R-L  circuits,  we  first  derived  a 
mechanical  analog  for  spring  and  damping  effects  from  [J.J.Blech,  “On  Isothermal  Squeeze 
Films,”  Journal  of  Lubrication  Technology,  v.  105,  1983,  pp.  615-620].  The  use  of  mechanical 
equivalent  circuit  models  had  previously  been  demonstrated  to  be  less  confusing  in  a  behavioral 
modeling  system  like  NODAS  that  has  both  electrical  and  mechanical  models. 
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The  model  contains  parallel  branches  of  series-connected  damper  and  spring  elements  (Figure 
35).  The  expressions  for  the  damper  and  spring  elements  are  where  m,  n  are  odd  integers,  w  is 
the  width  and  /  the  length  of  the  squeeze  film,  g  is  the  nominal  gap  which  is  equal  to  the 
thickness  of  the  squeeze  film.  Pa  is  the  ambient  pressure,  rjeff  is  the  effective  viscosity  of  the  air, 
and  z  is  the  displacement  of  the  oscillating  top  plate,  which  squeezes  the  air  film. 
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Figure  35.  (a)  Squeeze  film  between  two  parallel  plates,  (b)  Spring-damper  behavioral 

model. 

The  B-k  values  are  related  to  the  component  values  of  the  R-L  equivalent  circuit  model 
presented  in  [T.  Veijola,  H.  Kuisma,  J.  Lahdenpera,  and  T.  Ryhene,  “Equivalent-Circuit  Model 
of  Squeezed  Gas  Film  in  a  Silicon  Accelerometer,”  Sensors  and  Actuators  A,  vol.  48,  1995, 
pp.239-248]  as  B,„„  =  1/R,„n  and  ^,„„=  l/L^n.  The  velocity  (across  variable)  of  the  top  plate  was 
mapped  to  the  voltage  across  the  R-L  branches  and  the  damping  force  (through  variable)  was 
mapped  to  the  current. 

Once  the  existing  state  of  the  art  model  had  been  incorporated  into  NOD  AS,  we  began  a  series  of 
simulations  comparing  the  model  with  CFD-ACE+.  These  simulations  indicated  that  at  the  scale 
of  the  squeeze-film  found  in  the  differential  comb  finger  configurations  of  lateral  accelerometers 
(where  comb  fingers  are  2-5  p.m  thick),  the  behavioral  squeeze  film  model  was  off  by  as  much  as 
50%.  This  error  occurs  because  the  behavioral  model  inherits  Blech's  original  assumption  of 
trivial  boundary  conditions  that  set  the  gauge  pressure  to  zero  at  the  plate  edges.  This  assumption 
is  not  true  in  practice.  As  shown  from  the  pressure  distribution  obtained  by  CED-ACE-t-,  the 
gauge  pressure  is  zero  only  at  some  distance  from  the  plate  edges  (Figure  36). 
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Figure  36.  Pressure  distribution  on  the  surface  of  an  oscillating  plate  (Quarter  plate 

shown). 
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One  approach  to  handling  the  edge  effects  is  to  increase  the  plate  size  in  the  behavioral  model  by 
5L  to  match  the  numerical  simulation  results  with  non-trivial  boundary  conditions.  With  the 
extension  of  the  plate  dimensions,  the  effective  width  and  the  length  become  Wgji  =  w+^(w)  and 
=  l+dL(l).  These  effective  dimensions  are  used  in  equations  (4.1)  and  (4.2).  SL  depends  on  the 
geometry  of  the  squeeze  film  and  the  oscillation  frequency.  To  screen  the  importance  of  these 
factors,  3-D  numerical  simulations  with  frequencies  in  the  ranges  of  100  Hz  to  10  kHz,  gap  sizes 
from  1.5  )a.m  to  4  pm  and  plate  sizes  in  the  range  of  10  pm  to  1  mm  were  performed  and  the 
values  of  plate  extensions  ((5L)  needed  to  match  the  behavioral  and  numerical  results  determined. 
These  variable-screening  experiments  indicated  a  strong  dependence  of  SL  on  the  gap  and  the 
plate  size.  Simulations  were  performed  to  determine  the  dependence  of  &  on  these  variables. 
The  values  of  SL/g  that  gave  the  closest  match  between  the  behavioral  and  the  numerical 
simulations  are  determined  for  plate  sizes  in  the  range  of  10  pm  to  100  pm  (Figure  37).  Through 
a  linear  fit  minimizing  the  sum  of  squared  errors,  a  second-order  model  for  equivalent  plate 
extension  is  &(d)=g(0.8792+0.01d)  where  rf  is  a  variable  representing  /  or  w  in  microns.  The 
relation  indicates  that  SL  is  increasing  with  the  plate  size.  This  is  as  expected  because  larger 
plates  squeeze  more  air  and  cause  greater  pressure  perturbation  at  the  plate  edges.  The  pressure 
settles  to  ambient  pressure  further  away  from  the  plate  edges. 

Simulation  results  from  behavioral  models  with  and  without  inclusion  of  plate  extension  for  plate 
sizes  varying  from  10  pm  to  100  pm  are  compared  to  the  FEA  results.  Figure  38  plots  the  error 
in  the  two  behavioral  models.  The  model  with  plate  extensions  is  always  more  accurate.  The 
error  of  the  original  model  increases  drastically  for  plate  sizes  less  than  60  pm  on  a  side,  and  is 
still  about  10%  for  large  plate  sizes  due  to  the  trivial  boundary  conditions  assumed  in  its 
derivation.  The  enhanced  plate-extended  model  matches  the  FEA  results  to  within  3%  for  all  the 
plate  sizes  in  the  range  of  10  pm  to  100  pm. 

These  results  were  summarized  in  [S.  Vemuri,  G.  K.  Fedder,  T.  Mukherjee,  "Low  order  squeeze 
film  model  for  simulation  of  MEMS  devices",  Proc.  MSM,  2000.]  and  form  the  current  NODAS 
squeeze  film  model. 


Figure  37.  5L/gap  values  that  best  match  the  behavioral  and  numerical  results  for 

different  plate  sizes. 
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Accuracy  comparison  of  the  models 
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Figure  38.  Accuracy  of  the  modified  model  as  compared  to  the  original  one  for  different 

plate  sizes. 

The  final  verification  step  involved  experimental  verification  using  a  micromechanical  bandpass 
filter  fabricated  in  the  CMOS-MEMS  process  [Q.  Jing,  H.  Luo,  T.  Mukherjee,  L.  R.  Carley,  and 
G.  K.  Fedder,  “CMOS  Micromechanical  Bandpass  Filter  Design  Using  a  Hierarchical  MEMS 
Circuit  Library,”  Proc.  of  12th  IEEE  Inti.  Conf.  on  Micro  Electro  Mechanical  Systems  (MEMS 
’00),  Miyazaki,  Japan,  pp.  187-192,  January  23-27,  2000]. 

Figure  39  shows  the  NOD  AS  schematic  of  the  bandpass  filter.  The  filter  consists  of  3  crab-leg 
resonators,  with  differential  comb  transduction  of  the  input  signal  from  the  electrical  domain  to 
the  mechanical  domain  (on  the  left  resonator),  as  well  as  back  from  the  mechanical  domain  to  the 
electrical  domain  (on  the  right  resonator). 


resonator  resonator  spring  resonator 

Figure  39.  NODAS  schematic  of  a  bandpass  filter 
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The  frequency  response  of  the  filter  using  the  original  and  extended  NODAS  models  and  the 
experimental  results  are  shown  in  Figure  40.  The  simulation  results  are  in  excellent  agreement 
with  the  experimental  measurements.  The  results  show  that  the  error  in  computing  Q  due  to  the 
damping  model  was  reduced  from  20%  to  2%. 


-  Experimental 

NODAS  simulation  with  edge  effects 
NODAS  simulation  without  edge  effects 

Figure  40.  Output  voltage  of  a  micromechanical  CMOS  bandpass  filter. 

4.1.3.2  Parameterized  Models  of  Lateral  Damping  in  MEMS 

As  with  the  squeeze-film  damping,  this  sub-task  also  involved  four  stages:  we  first  implemented 
an  existing  lateral-damping  model;  then  we  compared  the  model  with  CED-ACE+;  and  studied 
the  comparisons  to  determine  how  best  to  extend  the  model;  finally  we  verified  the  extended 
model  against  both  CFD-ACE-i-  and  experimental  measurements. 

When  a  planar  structure  oscillates  laterally  over  an  immobile  substrate  plane,  the  resulting  gas 
flow  exerts  a  damping  force.  The  damping  force  has  a  spring  component  that  is  in  phase  with  the 
displacement  and  a  damping  component  in  phase  with  the  velocity.  The  complex  damping 
coefficient  representing  the  force  is  given  by  the  equation  (4.9).  The  real  part  represents  the 
damping  force  and  the  imaginary  part  represents  the  spring  or  the  inertial  force  [T.  Veijola, 
“Compact  Damping  Models  for  Lateral  Structures  Including  Gas  Rarefaction  Effects”,  Proc. 
MSM ,  20] 


n,ftAdv(z) 
Vo  dz 


z=0 


^eff  ^9 


cosh(qg)  +  cjksinh(qg) 

(1  +  q^X^  )sinh(qg)  +  2c(kcosh(qg) 


(4.9) 


where  Tjeff  is  the  effective  viscosity  coefficient  of  the  gas,  A  is  the  plate  area,  g  is  the  gap  between 
the  plate  and  the  substrate,  A,  is  the  molecular  mean  free  path  of  gas  at  the  operating  pressure  and 


48 


temperature,  and  q  is  the  complex  frequency  given  by 
the  plate  oscillation  frequency. 


with  p  is  air  density  and  co  is 


The  frequency  dependent  transfer  function  is  implemented  in  Verilog-A  as  a  Laplace  transform 
[S.  Vemuri,  “Behavioral  Modeling  of  Viscous  Damping  in  MEMS,”  Master's  Thesis,  Carnegie 
Mellon  University,  August  2000].  The  one  dimensional  slide  film  model  has  two  pins  indicating 
the  displacement  of  the  oscillating  plate  and  the  substrate.  This  is  converted  into  relative  velocity 
using  time  derivative  operator.  This  relative  velocity  is  used  to  compute  the  damping  force  in 
frequency  domain  and  inverse  Laplace  transformed  back  to  damping  force  in  time  domain.  The 
input  parameters  include  the  length  and  the  width  of  the  plate,  the  gap  and  the  ambient  pressure. 


To  include  the  edge  effects  in  the  behavioral  model,  the  length  and  the  width  of  the  plate  are 
each  increased  by  4  pm.  To  study  the  accuracy  of  the  model,  we  compare  the  behavioral 
simulation  results  with  the  finite  element  simulation  results.  As  the  frequency  domain  model  is 
intended  for  use  in  non-linear  transient  analysis,  a  comparison  of  square  wave  displacement 
patterns  is  reported  below.  The  first  three  dominant  harmonics  of  a  square  wave  as  the  input 
displacement  for  a  100  pm  square  plate  and  a  2  pm  gap.  The  damping  force  resulting  from  the 
oscillating  square  wave  displacement  is  shown  to  be  in  good  agreement  with  CFD  simulations  in 
Figure  42.  The  bottom  of  the  slide  film  is  considered  to  be  the  immobile  substrate. 


Figure  42.  Slide-film  model  comparison  for  an  oscillating  square  wave  displacement  on  the 

top-plate. 
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To  study  the  behavioral  model  performance  with  varying  frequency,  the  magnitude  and  the  phase 
of  the  damping  force  with  respect  to  the  plate  velocity  was  obtained  for  an  100  |J,m  square  plate 
and  a  2  pm  gap.  Figure  43  shows  the  magnitude  of  the  damping  force  Figure  44  shows  the  phase 
of  the  damping  force,  both  indicating  that  the  lumped  extended  plate  model  matches  well  with 
continuum  simulation. 


Figure  43.  Magnitude  of  peak  force  with  respect  to  oscillating  frequency  of  the  top  plate. 


Figure  44.  Phase  of  damping  force  with  respect  to  oscillation  frequency  of  top  plate. 

For  experimental  comparison,  we  verified  the  model  against  the  lateral  folded-flexure  resonators 
synthesized  using  the  resonator  synthesis  module  [T.  Mukherjee,  S.  Iyer,  and  G.  K.  Fedder, 
“Optimization-based  synthesis  of  microresonators,”  Sensors  and  Actuators  A,  70  (1998),  pp 
118-127]  developed  as  part  of  the  DARPA  Composite  CAD  BAA  96-16  “Foundations  for 
MEMS  Synthesis”  project.  The  measured  quality  factors  indicated  an  error  of  20%  with  the 
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original  model  with  the  error  reducing  to  8%  with  the  extended  model.  Unlike  the  squeeze  film 
model  where  extensive  CFD-ACE+  simulations  were  used  to  first  screen  for  the  best  variables, 
in  the  lateral  damping  case  we  extrapolated  the  results  from  the  squeeze-film  model  to  determine 
fixed  plate  extension  length  and  width  parameters.  The  CFD-ACE+  suite  of  runs  needed  for 
parameter  fitting  could  not  be  completed  by  the  end  of  the  project,  so  a  fixed  plate  extension  of  4 
fim  was  used.  This  is  the  primary  reason  behind  the  lower  accuracy  of  the  lateral  damping 
model. 


4.2  High-Fidelity  and  Reduced  Models  of  the  Synthetic  .Tets 

The  synthetic  jets  were  fabricated  and  measured  at  GT,  and  CFDRC  in  this  project  was 
developing  models  of  the  devices,  providing  high-fidelity  simulation  results,  design  tools,  and 
reduced  model  generation  techniques.  The  research  included  development  of  reduced-order 
models  of  the  synthetic  jets:  single-cell  domain  model,  implemented  in  3-cell  and  4-cell  models 
of  the  entire  device,  and  their  application  in  large-scale  simulations  of  synthetic  jet  arrays,  for 
active  electronic  cooling  and  air  flow  over  an  aircraft  wing.  All  the  details  of  the  synthetic  jets 
fabrication  and  measurements  are  described  below. 

4.2.1  Synthetic  Jets  Fabrication  and  Measurements 

4.2.1.1  Research  Highlights 

•  Microfabricated  two  types  of  synthetic  jet  modulator  arrays 

-  Electrostatic  Microvalve-type  Modulators 

-  Electrostatic  Membrane-type  Modulators 

•  Utilized  New  Microfabrication  Technologies  for  Synthetic  Jet  Modulator  Arrays 

-  Dry  Etching  of  Silicon  Using  Bosch  Deep-Trench  Silicon  Etch  Process  and  Plasma 
Therm  ICP 

-  Lamination  of  Dielectric  and  Metal  Layers  onto  Silicon  Substrates 

-  Polymer  Sacrificial  Layer  Removal  by  Two-Sided  Barrel  Etch  Step 

•  Explored  New  Characterization  Techniques  for  Synthetic  Jet  Modulator  Arrays 

-  Particle  Image  Velocimetry 

-  Infrared  Imaging  of  Heated  Surface  Cooled  by  Synthetic  Jet  Device 

•  Compiled  MEMS  Materials  Database  for  Device  Modeling/Simulation 

4.2.1.2  Introduction  to  Synthetic  Jets 

Conventional  jets  are  formed  by  fluid  flowing  from  an  area  of  high  pressure  to  an  area  of  lower 
pressure.  In  the  example  in  Figure  45,  a  pressure  drop  exists  across  an  orifice  plate  resulting  in  a 
jet  flow  from  left  to  right,  as  drawn,  so  long  as  pressure  Pi  >  P2. 
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Figure  45.  Conventional  jet  flowing  through  orifice  plate  due  to  pressure  drop  6P  =  P1-P2. 

Unlike  a  conventional  jet,  a  synthetic  jet  is  formed  by  an  oscillatory  flow  that  is  acoustically 
induced  near  the  edge  of  an  orifice.  A  synthetic  jet  actuator  consists  of  a  fixed  actuator  cavity 
bound  on  one  side  by  a  flexible  membrane  and  on  the  other  by  an  orifice  (Figure  46).  When  the 
membrane  is  vibrated  rapidly,  air  is  repeatedly  drawn  into  the  cavity  through  the  orifice  (Figure 
46a)  and  then  ejected  out  of  the  cavity  through  the  same  orifice  (Figure  46b).  As  the  outgoing 
flow  passes  the  sharp  edges  of  the  orifice,  the  flow  separates  forming  a  series  of  vortex  rings, 
which  propagate  normally  away  from  the  orifice  plate.  The  vortices  are  formed  at  the  excitation 
frequency,  and  the  nominally  round,  turbulent  jet  is  synthesized  by  the  interaction  of  the  vortices 
downstream  from  the  orifice.  Since  the  vortices  are  generated  by  an  oscillatory  flow  near  the 
edge  of  the  orifice,  the  actuator  transfers  momentum  to  the  air  without  net  mass  injection  into  the 
overall  system. 


Radial 

Distance 


Figure  46.  Cross  section  of  typical  synthetic  jet  actuator  showing  entrainment  of  ambient 
fluid  2a  (Pcavity  <  PAmbient)  and  vortex  generation  2b  (Pcavity  >  PAmbient)*  Figure  2b  also  shows 
a  typical  streamwise  velocity  profile  as  measured  downstream  of  the  orifice  plate. 
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4.2.2  Modulation  of  Synthetic  Jets  Using  Micromachined  Modulator  Arrays 

Three  synthetic  jet  modulation  schemes  have  been  the  focus  of  our  investigations.  In  the  direct 
approach,  an  array  of  micromachined  membrane  actuators  is  used  to  generate  the  synthetic  jet 
flow  (Figure  47). 

A  Jet  Flow 


□  Silicon  Substrate 

□  Polymer  Membrane 
CD  Metal  Electrode 

Figure  47.  Direct  generation/modulation  of  synthetic  jets  using  micromachined  membrane 

actuators. 

Modulation  of  the  jet  flow  is  achieved  by  varying  the  electrostatic  drive  signal  sent  to  the 
membrane  actuators.  So,  across  an  array  of  such  devices,  drive  signal  amplitudes,  frequencies, 
and  phasing  can  be  varied  giving  the  most  flexibility  in  terms  of  tuning  synthetic  jet  array  output. 
A  disadvantage  of  this  approach  is  that  to  obtain  that  flexibility,  a  high  voltage  amplifier  and  an 
oscillator  circuit  are  required  for  each  element  of  the  microjet  array. 

Since  high  jet  velocities  can  be  difficult  to  achieve  with  a  micromachined  electrostatic  membrane 
actuator,  a  variation  of  this  method  uses  the  micromachined  membrane  actuators  as  modulators 
of  jet  flow  generated  by  another  source,  typically  a  piezoelectrically-driven  membrane  as  in 
Figure  48. 


Membrane 

Dielectric  Polymer  Modulator 


Figure  48.  Schematic  of  membrane  modulation  scheme  (not  drawn  to  scale). 
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In  this  configuration,  the  piezoelectric  membrane  pneumatically  drives  the  membranes  of  all 
microjet  array  elements  in  parallel,  so  only  one  high  voltage  amplifier  is  required.  To  modulate 
an  individual  array  element,  a  DC  voltage  can  be  applied  to  the  membrane  electrode  of  that 
microjet  to  electrostatically  increase  the  apparent  stiffness  of  that  membrane,  reducing  the 
overall  membrane  movement  for  a  given  pneumatic  pressure.  Clamping  of  the  membrane  across 
the  orifice  hole  is  not  required  to  completely  shut  off  jet  flow;  instead,  reducing  the  membrane 
deflection  below  the  range  of  motion  required  for  synthetic  jet  generation  is  sufficient.  The 
polyimide  membrane  also  reduces  the  chance  of  arcing  by  more  completely  isolating  the  actuator 
electrode  from  the  substrate. 

The  third  modulation  method  investigated  utilizes  an  array  of  micromachined  valves  to  modulate 
the  output  of  a  synthetic  jet  generated  by  a  piezoelectrically  driven  membrane  (Figure  49). 


Figure  49  Schematic  of  electrostatic  micromachined  valve  modulator  structure  (not  drawn 

to  scale). 

The  microvalve  modulators  function  as  normally-open  valves  which  close  upon  application  of  a 
sufficiently  high  potential  difference  between  the  actuator  and  the  substrate.  For  actuation  to 
occur  the  gap  between  the  actuator  and  the  substrate  must  be  relatively  small,  otherwise  the 
electrostatic  force  will  be  too  small  to  close  the  valve.  The  small  gap  distance  produces 
significant  flow  degradation  when  comparing  similar  size  orifice  holes  with  and  without 
modulator  valves  positioned  over  them.  In  addition  to  overcoming  the  inherent  stiffness  of  the 
microvalve  suspension,  the  applied  voltage  must  also  be  able  to  overcome  the  force  of  the  air 
flowing  in  and  out  of  the  orifice.  If  the  applied  voltage  becomes  too  large,  arcing  can  occur 
which  will  permanently  damage  the  valve.  Again,  a  high  voltage  amplifier  supplies  the  drive 
signal  for  synthetic  jet  generation.  A  DC  voltage  applied  between  the  metal  valve  structure  and 
the  substrate  ground  plane  deflects  the  valve  plate  towards  the  substrate,  closing  the  valve. 

4.2.3  Fabrication  of  Synthetic  Jet  Modulators 

Since  the  majority  of  the  experimental  work  has  focused  on  the  development  and  testing  of 
microvalve  modulators,  we  begin  by  discussing  the  fabrication  of  the  silicon  microvalve 


54 


modulators.  Three  methods  have  been  investigated  for  fabricating  the  silicon  microvalve 
modulator  arrays  (Figure  50).  The  traditional  method  utilizes  standard  spin  coating  and  sputter 
deposition  techniques  to  deposit  the  polyimide  and  metal  layers,  respectively.  The  advantages  of 
the  traditional  method  include  good  control  over  layer  thickness  as  well  as  the  use  of  standard 
cleanroom  processing  equipment.  The  lamination  method  offers  quick  assembly  and  the  ability 
to  obtain  thicker  layers  than  typically  achieved  using  conventional  approaches.  All  three 
methods  utilized  the  Bosch  deep  silicon  trench  etching  process  to  etch  orifice  holes  through  the 
silicon  substrate  leaving  near  vertical  sidewalls.  Using  oxygen  based  plasma,  the  samples  are 
barrel  etched  to  remove  polyimide  located  under  the  actuator  and  its  moving  supports.  Figures 
52-55  show  several  examples  of  actuators  fabricated  using  the  lamination  method. 

Not  pictured  in  Figure  50  (due  to  its  similarity  to  the  traditional  method)  is  the  electroplating 
technique  in  which  the  sputter  deposition  step  is  used  for  seed  layer  deposition.  A  negative 
image  of  the  actuator  photolithography  mask  is  used  to  create  a  mold  of  thick  photoresist,  and 
electroplating  is  used  to  fill  that  mold  with  an  appropriate  metal,  typically  a  Ni/Fe  alloy.  Prior  to 
the  release  etch,  the  mold  and  exposed  seed  layer  are  removed.  The  main  disadvantage  of  this 
method  was  the  occasional  delamination  of  the  electroplated  plated  metal  due  to  the  heat 
generated  during  the  barrel  etching  process. 
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Figure  50.  Microvalve  modulator  fabrication  using  both  traditional  and  lamination 

methods. 


Figure  51.  Photograph  of  one  microvalve  of  a  modulator  array  fabricated  using  traditional 

fabrication  techniques. 
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Figure  54.  Photograph  of  released  aluminum  laminated  microvalve  modulator  arrays 
showing  several  of  the  actuator  designs  being  investigated. 


Figure  55.  Photograph  of  released  microvalve  modulators  fabricated  using  laminated 

copper  films. 
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The  membrane  modulator  fabrication  procedure  is  based  upon  the  process  used  to  fabricate  the 
microvalve  modulators  (Figure  56).  These  modulators  are  fabricated  as  two  separate  wafers, 
which  are  subsequently  bonded  together  with  epoxy.  Separate  fabrication  permits  spin  coating 
of  the  membrane  prior  to  the  silicon  etch  or  lamination  of  the  polymer  membrane  after  the  silicon 
etch  step.  A  tight  seal  between  the  two  wafers  is  achieved  by  clamping  the  cavity  and  orifice 
chips  together  with  magnets  prior  to  application  of  epoxy.  Figure  57  shows  a  photograph  of  a 
membrane  modulator  array  fabricated  using  the  method  described  in  Figure  56. 
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Figure  56.  Fabrication  sequence  for  membrane  modulator  arrays. 
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(a)  (b) 

Figure  57.  Photographs  showing  the  back  (a)  and  front  (b)  of  a  membrane  modulator 
array  fabricated  using  the  method  described  in  Figure  55. 

It  is  important  to  note  that  in  all  three  fabrication  procedures,  features  on  the  front  surface  of  a 
wafer  must  be  aligned  to  features  on  the  back  of  the  same  wafer  and,  in  some  cases,  to  features 
on  a  second  wafer.  This  makes  the  use  of  alignment  markers  critical  during  photolithography 
and  assembly  steps. 

Before  the  modulator  array  samples  can  be  characterized,  they  must  be  packaged  so  that 
electrical  contact  can  be  made  without  compromising  the  airflow  characteristics  of  the  device. 
For  the  membrane  type  devices  designed  to  generate  the  synthetic  jet  directly,  a  combination  of 
epoxy  and  conductive  epoxy  are  used  to  secure  wires  to  the  substrate  and  make  electrical  contact 
between  the  wires  and  the  bonding  pads,  respectively  (Figure  57). 


Figure  58.  Photograph  of  a  2x2  aluminum-laminated  microvalve  modulator  array 
mounted  on  a  patterned  copper-clad  PC  board  for  deflection  and  flow  testing. 

For  the  microvalve  and  membrane  modulator  samples,  PC  boards  with  patterned  copper  traces 
were  prepared  as  sample  mounting  plates  (Figure  58).  The  silicon  samples  are  placed  over  a 
large  hole  drilled  through  the  board  and  bonded  with  standard  epoxy.  Electrical  connections  are 
then  made  using  conductive  epoxy  from  the  bond  pads  on  the  silicon  sample  to  the  patterned 
copper  traces  on  the  test  board  Off-board  connections  can  then  be  made  merely  be  soldering 
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wires  directly  to  the  copper  traces  on  the  board.  As  shown  in  Figure  58,  the  board-mounted 
sample  is  ready  for  deflection  testing. 


For  flow  testing  the  sample  mounting  boards  may  be  attached  to  the  front  of  the  piezoelectric 
driver  mounting  plate.  (Figure  59).  To  ensure  an  airtight  seal  between  the  piezo  mounting  plate 
and  the  sample  mounting  plate,  a  rubber  o-ring  is  used.  Figure  60  shows  a  schematic  of  the 
piezoelectric  membrane  actuator  used  to  generate  the  synthetic  jets  when  testing  the  microvalve 
and  membrane  modulator  samples.  The  driver  can  be  bonded  to  the  piezo  mounting  plate  by  a 
variety  of  adhesives. 
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Back  View  Of  Piezo  Mounting  Plate 


I  I  1/4"  Aluminum 


Cross  Section  After  Assembly 


Compressed  O-Ring  (DASH#220) 
Piezo  Ceramic  Diaphragm 
Silicon  Modulator  Array  Chip 


Top  View  of  Sample  Mounting  Plate 


I  I  1/16"  FR4  Copper  Clad  Circuit  Board 
Sample  Mounting  Area 


Figure  59.  Schematic  showing  details  of  the  piezo  driver  and  sample  mounting  plates. 
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Piezo  Ceramic  Diaphragm 


Murata  7NB-31R2-1  Specifications 

Resonant  Frequency  =  1.3kHz 
Resonant  Impedance  =  300  ohm 
Capacitance  (@1?0Hz)  =  4DnF 


Operating  Conditions  During  Testing 

Drive  Frequency  =  1470  Hz 
Drive  Voltage  =  40Vrms  sine  wave 


I  I  Silver  Electrode 


Piezo  Ceramic 


Figure  60.  Schematic  showing  details  of  the  piezoelectric  driver  element  used  in  these 

experiments. 

4.2.4  Characterization  of  Synthetic  Jet  Modulator  Arrays 


Prior  to  flow  testing  of  the  modulator  arrays,  deflection  testing  was  performed  to  determine  the 
minimum  voltages  required  to  close  the  microvalve  actuators  or  to  visibly  deflect  the  polymer 
membranes  of  the  membrane  modulators.  Mounted  samples  were  photographed  using  a  high- 
resolution  video  microscope  as  shown  in  Figure  61.  By  carefully  adjusting  the  angle  of  the 
incident  light  a  and  the  angle  of  the  sample  relative  to  the  camera  lens  P,  a  high-resolution,  high 
magnification  image  with  minimal  reflections  could  be  obtained. 


Figure  61  Deflection  testing  setup. 
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When  tested  without  the  synthetic  jet  actuator  attached,  microvalve  modulator  array  elements 
have  been  shown  to  actuate  using  applied  voltages  as  low  as  100  volts.  The  lowest  actuation 
voltages  measured  for  the  membrane  modulators  were  significantly  higher  at  500  volts. 

Once  the  minimum  actuation  voltage  was  determined,  one  or  more  flow  visualization  techniques 
were  used  to  determine  if  the  synthetic  jet  could  be  modulated  using  voltages  at  or  above  the 
minimum  actuation  voltage  previously  measured.  The  simplest  visualization  technique  used  is 
known  as  smoke  visualization.  Since  synthetic  jets  are  composed  entirely  of  entrained  ambient 
air,  a  smoke  source  placed  near  the  orifice  of  the  jet  will  allow  the  jet  to  entrain  and  eject  smoke 
making  the  synthetic  jet  flow  visible  (Figure  62).  Once  the  airflow  is  visible,  one  may  visibly 
determine  if  modulation  has  been  achieved. 


Two  types  of  smoke  were  used  in  these  experiments.  A  commercially  available  liquid 
commonly  used  in  theatrical  fog  machines  can  be  heated  electrically  until  smoke  is  produced. 
Smoldering  incense  was  also  used  to  generate  smoke  for  jet  entrainment.  Since  the  smoke  from 
either  of  these  heated  sources  will  rise,  the  sample  was  suspended  vertically  above  the  smoke 
source  so  that  a  thin  sheet  of  smoke  would  rise  past  the  sample’s  orifices.  In  practice,  both  types 
of  smoke  were  found  to  clog  the  synthetic  jet  orifices  and  modulators  quickly.  Wherever  the  fog 
machine  smoke  touched  the  sample,  it  left  an  oily  residue,  which  proved  to  be  difficult  to  remove 
without  destroying  the  sample.  This  residue  cut  the  testing  lifetime  to  a  mere  5-10  minutes  in 
some  cases.  The  smoke  particles  from  the  incense  also  clogged  the  microjets,  though  at  a  slower 
rate.  Another  side  effect  of  smoke  visualization  observed  in  earlier  experiments  was  that  the 
updraft  caused  by  the  rising  smoke  deflected  the  synthetic  jets  upwards  so  that  they  were  no 
longer  normal  to  the  surface  of  the  silicon  chip.  As  such  smoke  visualization  was  found  to  be  an 
inadequate  method  for  testing  the  synthetic  jet  modulator  arrays. 


Synthetic  ^ 
Jets  <■ 


Sample 

Mounting 

Plate 


Smoke 

Source 


Figure  62.  Smoke  visualization  of  synthetic  jet  actuators. 

A  Schlieren  optical  system  was  also  used  to  quickly  determine  if  the  modulator  arrays  were 
modulating  the  synthetic  jets.  As  shown  in  Figure  63,  a  Schlieren  optical  system  can  be  used  to 
visualize  small  index  of  refraction  changes  occurring  in  the  test  area.  The  light  passing  through 
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the  material  with  the  different  index  of  refraction  refracts  down  a  slightly  different  path  where  it 
is  blocked  by  a  knife  edge  producing  a  shadow  in  the  image  seen  by  the  video  camera. 


To  use  this  non-destructive  means  of  flow  visualization,  we  had  to  introduce  an  index  of 
refraction  change  in  the  synthetic  jet.  One  way  to  do  this  is  to  heat  the  air  of  the  jet  being  tested. 
Previous  experiments  had  shown,  however,  that  heating  of  the  jets  would  not  only  introduce 
thermal  drafts  that  could  obscure  the  synthetic  jets,  but  the  heat  could  also  interfere  with  the 
operation  of  the  microvalves. 

Another  way  to  make  the  synthetic  jet  visible  is  to  allow  it  to  entrain  a  gas  with  an  index  of 
refraction  different  than  that  of  air.  Helium  entrainment  produces  an  image  with  reasonable 
contrast,  but  we  determined  experimentally  that  the  breakdown  voltage  of  the  electrostatic 
actuator  was  significantly  lower  in  a  helium  rich  atmosphere  than  in  the  normal  ambient  air 
environment  used  for  measuring  the  deflection  voltage.  So,  many  microvalves,  which  would 
actuate  in  a  normal  air  environment,  would  experience  catastrophic  failure  during  the  Schlieren 
visualization  experiment.  In  an  attempt  to  solve  this  problem,  carbon  dioxide  gas  was  also  used. 
Since  the  breakdown  voltage  of  a  carbon  dioxide-rich  atmosphere  is  higher  than  that  of  normal 
air,  modulators  which  deflect  in  normal  air  without  breakdown  will  also  deflect  in  the  carbon 
dioxide-rich  atmosphere  without  breakdown.  The  index  of  refraction  change,  however,  is  so 
small  that  the  contrast  in  the  images  between  the  ambient  air  and  the  synthesized  jet  was  too 
small  most  times  to  determine  if  the  modulator  arrays  were  operational.  We  also  observed  that 
the  entrainment  of  the  visualization  gases  could  distort  the  synthetic  jet  as  it  emerged  from  the 
sample  making  the  jet  non-perpendicular  to  the  surface  of  the  silicon  chip. 

An  additional  difficulty  in  applying  the  Schlieren  visualization  technique  is  that  the  image 
contrast  of  the  synthetic  jets  is  reduced  by  the  small  size  of  the  jets  we  were  investigating.  When 
imaging  a  larger  jet,  say  a  centimeter  scale  or  larger  jet,  the  overall  distance  that  the  light  must 
travel  within  the  material  of  different  refractive  index  is  physically  longer  than  the  path  traversed 
through  the  sub-millimeter  scale  synthetic  jet.  The  longer  path  generates  more  beam 
displacement  at  the  knife  edge  and,  hence,  more  contrast  in  the  image.  Another  disadvantage  of 
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using  this  technique  to  visualize  the  operation  of  the  modulator  arrays  is  that  the  only 
information  this  technique  can  provide  is  quick  visual  confirmation  of  modulator  operation.  No 
velocity  information  or  any  other  data  is  gathered  by  this  technique. 

Other  non-destructive  methods  have  also  been  used  to  test  the  modulator  arrays  with  various 
degrees  of  success.  Pitot  tube  velocity  measurements  have  been  made  of  the  microjet  devices. 
For  these  measurements,  a  small  diameter  tube  is  connected  to  a  baritron  by  thin  tubing  (Figure 
64).  The  tube  is  positioned  in  front  of  the  jet  orifice  such  that  the  jet  blows  into  the  tube.  For 
jets  with  relatively  high  velocities,  the  measurements  of  streamwise  velocity  are  accurate 
provided  that  the  pitot  tube  itself  is  far  enough  from  the  orifice  that  the  size  of  the  tube  is  not  of 
the  same  order  of  magnitude  as  the  size  of  the  jet  being  measured.  For  low  velocity  jets,  the 
accuracy  of  the  measurement  is  limited  by  the  sensitivity  of  the  baritron  and  by  the  conductance 
of  the  tubing. 


Figure  64.  Pitot  tube  velocity  measurement  (A)  good  pitot  placement  and  (B)  pitot  too  close 

for  accurate  measurement. 


A  hot-wire  anemometer  has  also  been  used  for  non-destructive  testing  of  the  synthetic  jet 
modulators.  In  these  experiments,  the  sample  is  mounted  to  a  fixed  stage  and  an  x-wire  probe  is 
attached  to  an  x-y-z  positioning  stage.  Under  computer  control,  the  probe  may  be  scanned  across 
the  sample  to  record  2-axis  velocity  data.  The  main  disadvantages  experienced  using  this 
technique  were  related  to  calibration  and  probe  size.  As  built,  our  test  setup  was  geared  to 
testing  larger,  higher  velocity  jets,  and  the  smaller  velocities  of  the  microvalve  modulator  jets 
were  at  or  below  the  measurable  range  of  the  anemometer.  As  in  the  case  of  the  pitot  tube,  the 
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probe  size  was  also  on  the  order  of  the  size  of  the  jet  being  tested  making  the  accuracy  of  the 
measurements  questionable,  especially  near  the  sample  surface  where  the  highest  velocities  were 
found.  Finally,  the  x-wire  probe  was  an  impractical  choice  for  our  final  jet  impingement 
experiments  (described  later  in  this  section)  since  the  probe  would  block  placement  of  the 
impingement  plate. 

The  best  technique  we  found  for  gauging  the  performance  of  the  micromachined  synthetic  jet 
modulators  is  known  as  particle  image  velocimetry  (PIV),  To  use  this  technique,  a  cross  section 
of  the  flow  field  was  illuminated  by  a  thin  sheet  of  laser  light  (Figure  65).  A  high  speed,  high 
resolution  video  camera  was  then  used  to  photograph  particles  suspended  in  the  flow  field. 
Given  the  length  scale  of  the  flow  and  the  sampling  time,  image  processing  software  can 
determine  the  velocity  vectors  associated  with  various  points  in  the  image.  As  in  the  case  of 
smoke  visualization,  the  smoke  particles  used  may  clog  the  jets,  but  the  PIV  technique  obtains  a 
complete  velocity  map  for  the  points  illuminated  by  the  laser  sheet.  Moreover,  this  technique 
allows  velocity  data  to  be  recorded  during  the  jet  impingement  experiment. 

Top  View  of  PIV  Test  Setup 
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Currently  under  investigation  is  a  new  imaging  technique  in  which  jet  modulation  is  gauged  by 
observing  the  cooling  effect  of  the  modulated  synthetic  jet  upon  a  heated  object.  By  examining 
infrared  images  of  the  cooled  object,  we  hope  to  be  able  to  demonstrate  three-dimensional 
operation  of  the  modulator  array.  As  shown  in  Figure  66,  the  copper-clad  Kapton  substrate  is 
heated  by  a  thin-film  Kapton  heating  element  attached  to  the  copper  side  of  the  substrate.  The 
synthetic  jet  is  positioned  such  that  it  will  impinge  upon  the  heating  element.  An  infrared 
camera  is  used  to  view  the  temperature  distribution  of  the  substrate  from  the  Kapton  side.  A 
major  advantage  of  this  test  method  is  that  it  is  non-destructive  in  nature,  requiring  no  smoke  or 
visualization  gas  entrainment  that  could  interfere  with  the  operation  of  the  modulator  array. 

To  calibrate  the  temperature  measurement  setup,  an  accurate  emissivity  value  is  needed  for  the 
copper-clad  Kapton  film.  To  determine  the  unknown  emissivity,  a  material  of  known  emissivity 
is  used.  A  small  piece  of  thin  black  tape  with  an  emissivity  of  0.95  was  placed  in  the  center  of 
the  heated  substrate  facing  the  infrared  camera.  With  a  constant  power  of  0.9W  delivered  to  the 
heating  element,  the  temperature  of  the  black  tape  was  allowed  to  stabilize,  and  the  temperature 
of  the  tape  was  measured  by  the  infrared  camera  using  the  known  emissivity  of  0.95  for  the  tape. 
The  tape  was  then  removed,  and  the  temperature  of  the  heated  substrate  was  allowed  to  stabilize 
once  again.  The  emissivity  was  then  recalibrated  to  make  the  measured  temperature  of  the 
substrate  agree  with  the  measured  temperature  of  the  black  tape. 
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Infrared  Testing  Setup 
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Figure  66.  Schematic  of  infrared  test  setup  used  to  evaluate  the  cooling  effect  of  the 
modulated  synthetic  jet  arrays  on  a  heated  surface. 

4.2.5  Experimental  Data 

Given  the  volume  of  data  collected,  we  first  include  a  sample  of  the  data  collected  during  the 
characterization  of  a  typical  microvalve  modulator  array.  The  data  presented  is  from  a  single 
microvalve  modulator  array  (as  pictured  in  the  top  left  comer  of  (Figure  55)  mounted  for  testing 
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as  shown  in  Figure  59  with  the  piezoelectric  driver  element  picture  in  Figure  60  used  for 
synthetic  jet  generation.  Initial  deflection  testing  was  performed  using  the  video  microscope 
setup  pictured  in  Figure  61  with  modulation  voltages  from  406V  to  707V  applied.  Initial  flow 
testing  was  performed  using  the  pitot  tube  setup  diagramed  in  Figure  64.  Jet  velocities  from 
7m/sec  for  all  four  jets  operating  to  9.7m/sec  for  a  single  jet  operating  were  recorded  using  the 
pitot  tube  measurement  apparatus  with  a  42Vrms,  1470Hz  sine  wave  drive  signal  supplied  to  the 
piezoelectric  driver.  Subsequent  testing  was  performed  with  a  40Vrms,  1470Hz  sine  wave  to 
extend  the  lifetime  of  the  piezoelectric  driver.  Temperature  distributions  of  the  synthetic  jets 
impinging  on  a  heated  surface  were  measured  using  the  infrared  camera  setup  sketched  in  Figure 
66,  and  corresponding  velocity  data  for  this  experiment  were  collected  using  the  PIV 
measurement  setup  pictured  in  Figure  65. 


Actuator  Layout  Diagram 
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Figure  67.  Schematic  showing  mask  patterns  and  layer  thickness  of  typical  synthetic  jet 

microvalve  modulator  array. 

Figures  68,  69,  70  and  71  show  vorticity  plots  for  the  two  jets  in  the  flow  field  illuminated  by  the 
PIV  laser  sheet  (Figure  65).  In  these  Figures,  the  synthetic  jet  modulator  array  flow  is  impinging 
upon  the  heated  surface  as  pictured  in  Figure  66.  Figures  72, 73,  74  and  75show  PIV  streamwise 
velocity  (V)  profiles  for  the  Left  Jet,  Right  Jet,  Two  Jets,  and  Four  Jets  cases,  respectively. 
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Figures  76,  77,  78  and  79  are  PIV  velocity  profiles  showing  the  perpendicular  velocity 
component  U.  Note  that  the  sign  of  U  changes  depending  upon  the  direction  of  entrainment. 
Figures  80,  81,  82  and  83  are  PIV  speed  contours  of  the  four  cases  examined,  and  Figures  84,  85, 
86,  87  and  88  show  the  corresponding  temperature  distributions  for  no  cooling.  Left  Jet,  Right 
Jet,  Both  Jet,  and  Four  Jets,  respectively. 


Figure  68.  PIV  vorticity  plot  for  left  synthetic  jet  (called  “Left  Jet”)  from  microvalve 

modulator  array. 
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Figure  69.  PIV  vorticity  plot  for  right  synthetic  jet  (called  “Right  Jet”)  from  microvalve 

modulator  array. 


Figure  70.  PIV  vorticity  plot  for  two  synthetic  jets  (those  two  illuminated  by  the  laser  sheet 
known  as  “Both  Jets”)  from  microvalve  modulator  array. 
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Figure  71.  PIV  vorticity  plot  for  two  synthetic  jets  from  microvalve  modulator  array  with 
all  four  of  the  microvalves  open  (called  “Four  Jets”). 


Figure  72.  PIV  streamwise  velocity  profileV  for  just  the  left  synthetic  jet  from  microvalve 

modulator  array  (“Left  Jet”). 


Figure  73.  PIV  streamwise  velocity  profile  V  for  the  right  synthetic  jet  from  microvalve 

modulator  array  (“Right  Jet”). 


Figure  74.  PIV  streamwise  velocity  profile  V  for  the  two  synthetic  Jets  from  microvalve 
modulator  array  illuminated  by  the  laser  sheet  (“Both  Jets”). 
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Figure  75.  PIV  streamwise  velocity  profile  V  for  the  two  synthetic  jets  from  microvalve 
modulator  array  when  all  four  valves  are  open  (“Four  Jets”)* 
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Figure  76.  PIV  perpendicular  velocity  profileV  for  just  the  left  synthetic  jet  from 
microvalve  modulator  array  (“Left  Jet”). 
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Figure  77.  PIV  perpendicular  velocity  profileV  for  just  the  right  synthetic  Jet  from 
microvalve  modulator  array  (“Right  Jet). 


Figure  78.  PIV  perpendicular  velocity  profileV  for  both  synthetic  jets  from  microvalve 

modulator  array  (“Both  Jets”). 


Figure  80.  PIV  speed  contours  for  Left  Jet  case. 
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Figure  81.  PIV  speed  contours  for  Right  Jet  case. 


Figure  82.  PIV  speed  contours  for  Both  Jets  case. 
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Figure  83.  PIV  speed  contours  for  Four  Jets  case. 


Figure  84.  Inframetrics  ThermaCAM  image  of  temperature  distribution  of  heated  surface 

without  jet  cooling. 
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Figure  85.  Inframetrics  ThermaCAM  image  of  temperature  distribution  of  heated  surface 

as  cooled  by  the  left  jet  only. 


Figure  86.  Inframetrics  ThermaCAM  image  of  temperature  distribution  of  heated  surface 

as  cooled  by  the  right  jet  only. 
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Figure  87.  Inframetrics  ThermaCAM  image  of  temperature  distribution  of  heated  surface 

as  cooled  by  both  jets. 
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Figure  88.  Inframetrics  ThermaCAM  image  of  temperature  distribution  of  heated  surface 

as  cooled  by  all  four  jets. 

We  also  include  a  sample  of  the  data  collected  from  experiments  with  the  membrane  modulator 
devices. 

In  this  case  the  cavity  side  of  the  membrane  modulator  array  tested  is  pictured  in  Figure  89a,  and 
the  orifice  side  of  the  modulator  array  is  pictured  in  Figure  89b.  Figure  89c  shows  the 
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undeflected  membrane,  and  Figure  89d  shows  a  membrane  deflected  by  an  800V  actuation 
signal.  Figure  90  shows  the  results  of  the  jet  cooling  experiment  performed  with  this  device.  In 
this  test  a  lower  resolution  infrared  camera  (Hughes  Probeye)  was  used  to  visualize  the  heated 
surface.  From  Figure  90a  to  Figure  90b,  one  can  see  a  decrease  in  the  diameter  of  the  hotspot  as 
a  result  of  the  cooling  action  provided  by  the  membrane  modulator  device.  As  typical  for  the 
membrane  modulator  devices,  this  sample  later  failed  because  the  high  actuation  voltage 
required  shorted  out  the  actuators  before  PIV  velocity  data  could  be  obtained. 


(c)  (d) 

Figure  89.  Photographs  showing  actuation  of  membrane  modulator.  Assembled  sample 
Membrane  side  and  (b)  Orifice  side.(c)  Undefiected.  (d)  Deflected. 

Again  note  shadow  and  focal  changes. 
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(a)  (b) 

Figure  90.  Infrared  photographs  showing  the  cooling  of  a  heated  surface  with  the  synthetic 
jet  produced  by  a  membrane  modulator  device,  (a)  no  synthetic  jet.  (b)  synthetic  jet 
activated.  Note  the  decrease  in  diameter  of  the  center  hot  spot  upon  activation  of  the 

synthetic  jet. 

4.2.6  High-Fidelity  and  Reduced  Models 

Synthetic  jets  are  generated  at  the  orifice  of  a  partially  open  cavity  with  an  oscillating  membrane 
wall  opposite  to  the  orifice  opening  as  shown  schematically  in  Figure  91  [Smith  and  Glezer, 
1998].  Microfabricated  arrays  of  synthetic  jets  are  being  explored  for  applications  in 
aerodynamic  flow  control,  in  cooling  of  electronics  packages,  and  in  mixing  in  microchemical 
reactors. 


4  i 

Driver 

Figure  91.  Schematic  of  synthetic  jet. 

Figure  92.  High-fidelity  model  of  synthetic  jet. 

The  unsteady  flow  dynamics  inside  the  cavity  of  the  synthetic  jet  is  analyzed  with  CFD-ACE+. 
The  code  provides  multi-disciplinary  simulation  capability  by  coupling:  fluid  flow,  heat  transfer, 
mixing  and  chemistry,  stress/deformation,  electrofluidics,  electrostatics,  electromagnetics,  and 
other  discipline  field  solvers  specifically  adapted  for  MEMS  [Przekwas,  1999].  The  fluid  flow 
model  solves  the  time  dependent  continuity  equation,  the  pressure  based  Navier-Stokes  equations 
and  energy  balance  equation.  In  the  present  formulation  they  are  written  in  a  strong  conservation 
integral  form  on  time  dependent  arbitrary  moving/deforming  geometry  and  are  equally 
applicable  to  incompressible  and  compressible  flows.  Numerical  solution  of  this  equation  set 


requires  discretization  of  the  computational  domain  into  a  large  number  of  generalized  control 
volumes.  To  allow  full  freedom  in  the  control  volume  shapes,  the  governing  equations  are 
expressed  in  the  integral  flux  form: 


—  Jpdv  +'|p(v  -  Vg )  •  nda  =  0 


at 


d 
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3  3 
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where  Ui  is  the  ith  Cartesian  component  of  the  velocity,  Vg  is  the  grid  velocity  due  to  grid 
motion,  p  is  the  static  pressure,  ht  is  the  total  enthalpy,  Xy  is  the  stress  tensor  for  both  laminar  and 
turbulent  flows  and  fj  is  the  ith  Cartesian  component  of  body  force.  For  turbulent  flows,  the 
Reynolds  stress  tensor  is  closed  with  a  standard  k-E  model.  An  implicit  second  order  time-space 
control  volume  discretization  method  is  used  to  solve  the  above  equations  on  unstructured 
meshes. 

The  equations  are  solved  using  CFD-ACE+  software  developed  at  CFDRC  with  specific 
application  for  MEMS.  A  unique  feature  of  the  flow  solver  is  the  application  of  generalized 
polyhedra  control  volumes  for  super  resolution  of  jet  cavities  with  a  single  control  volume  with 
large  number  of  faces  (cavity  walls,  membrane,  orifice,  ...).  A  moving  wall  and  deforming  grid 
method  is  used  to  model  the  interior  of  the  cavity  and  a  free  space  environment  is  used  to  model 
the  jet  evolution  (Figure  92).  The  oscillating  membrane  is  periodically  ingesting  and  expelling 
the  air  from  the  cavity.  Several  oscillation  cycles  are  performed  to  obtain  cycle  independent 
results.  The  model  is  validated  against  the  experimental  data  obtained  GT  for  the  jet  actuated 
with  1  kHz  oscillation  frequency.  Unsteady  flow  field  have  been  experimentally  measured  using 
Particle  Image  Velocimetry  (PIV)  technique. 

The  height  of  the  actuation  cavity  is  95  mm,  cavity  width  is  19  mm  and  the  orifice  diameter  is 
6.35  mm.  Transient  computational  simulations  have  been  performed  with  membrane 
deformation  dynamics  represented  as  a  moving  wall  with  prescribed  harmonic  motion.  As  the 
membrane  dynamics  could  not  be  measured  in  the  experiments,  a  direct  comparison  between 
computational  and  experimental  data  was  not  possible.  Instead,  several  computational  runs  were 
performed  to  match  the  maximum  flow  rate  through  the  orifice  with  a  membrane  displacement  as 
a  parameter.  A  comparison  between  experimental  and  computational  vortical  flow  structures  at 
several  time  instances  is  presented  in  Figure  93.  Time  average  velocity  profiles  at  several  axial 
locations,  measured  at  GT  and  calculated  with  CFD-ACE  are  shown  in  Figure  94.  Considering 
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the  uncertainty  of  boundary  conditions  and  assumed  membrane  dynamics,  good  agreement 
between  predicted  and  experimental  velocity  profiles  has  been  achieved  (Figure  94). 


PIV  Measurements 


Figure  93.  Synthetic  jet  simulation  results  compared  to  measurements. 
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Figure  94.  Mean  velocity  profiles  at  different  height  position  over  the  synthetic  jet  orifice: 
a)  measured  with  PIV  technique  at  Georgia  Tech,  h)  obtained  from  CFD-ACE+  simulations  at 

CFDRC. 
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Reduced  Models  of  Synthetic  Jets 

Computational  modeling  of  large  numbers  of  synthetic  jet  cavities  with  high  fidelity  3-D  models 
is  not  practical  because  of  large  computational  cost.  A  "compact  model"  accounting  for  all 
essential  physics  expressed  in  terms  of  algebraic  equations,  or  ordinary  differential  equations,  is 
needed  to  model  large  arrays  of  jets.  There  are  several  approaches  to  formulate  a  compact  model 
including  analytical  solutions,  curve  fits  to  experiment^  data,  or  equivalent  electrical  circuits. 
Equivalent  circuit  models  of  air  dumping  for  inertial  sensors  have  been  recently  demonstrated  by 
the  authors  earlier. 

Single-Cell  Model  of  Synthetic  Jet 

In  this  project,  we  have  developed  a  novel  concept  of  mixed-dimensionality  approach  combined 
with  the  polyhedral  grid  capability  to  model  a  synthetic  jet  with  a  single  control  volume  or  a  one¬ 
dimensional  approximation.  Both  capabilities  have  been  implemented  in  CFD-ACE+  tool.  In 
the  proposed  approach,  to  represent  the  cavity,  a  single  control  volume  (CV)  fully  conforming  to 
the  cavity  geometry  is  used.  This  single  polyhedral  CV,  however,  has  a  large  number  of  cell 
faces,  some  of  them  are  fixed  wall  segments,  the  vibrating  membrane,  or  the  orifice  opening. 
Full  set  of  general  conservation  equations  (mass,  momentum,  energy)  is  solved  in  a  single  cell  in 
time-accurate  manner.  The  cavity  volume  is  linked  through  the  orifice  with  the  external  flow 
field  in  a  fully  implicit  manner  resulting  in  a  robust  simulation  algorithm.  If  a  desired  accuracy 
of  the  compact  model  is  achieved,  a  large  number  of  cavities  can  be  modeled  in  a  very  cost 
effective  manner. 

Figure  95  shows  the  geometry  of  a  reduced  model  of  the  synthetic  jet,  composed  of  three 
cylindrical  domains,  each  of  them  modeled  with  a  single  discretization  cell  having  32  boundary 
faces. 


Figure  95.  Single-cell  model  of  synthetic  jet  (3  domains  =  3  cells);  (a)  geometry;  (b)  velocity 
distribution,  calculated  using  the  3-ceIl  model  of  synthetic  jet. 
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Discussion  of  Reduced  Model  Validity  and  Accuracy 

Our  initial  studies  of  compact  models  of  synthetic  jets  showed  that  for  the  present  configuration 
a  two-cell  model  (chamber  and  orifice)  was  needed  for  the  cavity,  and  an  additional  domain  for 
the  external  environment  (surrounding  air).  Table  4.1  presents  relative  comparison  of 
computational  mesh  requirements  and  CPU  times  for  3-D,  ID  and  one-cell  models  to  simulate 
the  single  cycle  of  the  oscillatory  jet.  CPU  savings  for  a  single  cell  are  already  apparent  when 
compared  with  the  3-D  (full  cell)  model.  One  of  the  model  requirements  is  good  accuracy  for 
widest  possible  operating  range  (frequency,  displacement,  pressure  level,  geometry,  etc.). 
Computational  studies  showed  that  for  small  frequencies  (up  to  200  Hz)  incompressible  and 
compressible  flow  models  are  very  similar.  For  higher  frequencies,  large  discrepancies  in  flow 
rate  have  been  observed.  Figure  96  presents  both  characteristics  for  up  to  1  KHz.  We 
recommend  that  a  compressible  version  of  the  model  should  be  used  for  entire  range  of 
frequencies. 


Table  4.1.  Computational  mesh  resolution  and  CPU  Time  for  3-D,  ID  and  one  cell  models 
_ _ of  a  synthetic  jet. _ _ 


Model  Type 

Number  of  Discretization  Elements  (axial  x  radial)  in: 

CPU  Time 
(200  time 
steps) 

Chamber 

Orifice 

Outside 

Total 

Full  Cell 

64x34 

lOx  14 

64x38 

4748 

41  min  44  sec 

ID  Cell 

64  X  1 

10  X  1 

10  X  1 

84 

1  min  48  sec 

One  Cell 

_ iLl _ 

1  X  1 

1  X  1 

3 

23  sec 

At  low  frequencies,  the  pressure  field  inside  the  cavity  is  quite  uniform.  Above  300  Hz  we 
found  that  acoustic  (pressure  wave)  effects  became  important  and  a  single  cell  cavity  model 
would  not  capture  the  spatial  variation  of  pressure  along  the  axial  direction  within  the  cavity. 
For  high  frequency  applications,  we  have  proposed  a  one  dimensional  model  with  exactly  the 
same  equations  as  in  the  one-cell  model  but  solved  on  axially  arranged  polyhedral  elements. 


Maximum  Mass  Flux  at  Orifice 
Compressible  and  Incompressible  Flows 


Figure  96.  Results  of  3-D  high-fidelity  computations  of  mass  flow  in  synthetic  jet  orifice,  using 

compressible  and  incompressible  gas  models. 
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Figure  97  presents  the  comparison  of  predicted  mass  flow  rate  through  the  orifice  for  3-D,  ID 
and  one-cell  models  for  100  Hz  frequency.  Good  agreement  between  all  three  models  has  been 
observed. 


Mass  Flux  at  Orifice  (Compressible  Flow) 
Frequency  =  100  Hz 


Time  (sec,) 

Figure  97.  Mass  flow  calculated  using  various  fldelity  models  for  actuation  frequency  100  Hz. 


For  high  frequency  (1  KHz)  operation,  as  shown  in  Figure  98,  a  very  large  discrepancy  exists 
between  the  one-cell  model  and  ID,  3-D  models.  At  this  frequency,  pressure  wave  effects 
become  dominant,  and  a  ID  reduced  model  has  to  be  used.  Good  accuracy  agreement  between 
3-D  and  ID  models  has  been  achieved  without  any  "tuning"  of  the  ID  model. 


Mass  Flux  at  Orifice  (Compressibel  Flow) 
Frequency  =  1000  Hz 


Figure  98.  Mass  flow  calculated  using  various  fldelity  models  for  actuation  frequency  1000  Hz. 


87 


Conclusions 

We  have  developed  a  novel  concept  of  compact  models  for  synthetic  jets.  The  model  uses  a 
polyhedral  control  volume  capability  of  CFD-ACE+  software  tool  to  model  complex  dynamic  3- 
D  shapes  with  moving  walls  with  multiple  inlets/outlets  with  a  single  cell  "super  element".  The 
model  has  been  validated  against  3-D  high-fidelity  simulation  data  obtained  for  a  range  of 
parameters  such  as  geometry,  actuation  frequency,  amplitude,  operational  pressure,  etc. 

It  was  found  that  the  incompressible  form  of  the  model  was  valid  for  up  to  200  Hz  actuation. 
Large  errors  were  observed  at  higher  frequencies.  Compressible  formulation  is  more  general  and 
accurate  for  full  range  of  frequencies.  At  high  frequencies,  where  pressure  wave  effects  become 
important  within  the  cavity,  one  cell  model  has  to  be  replaced  by  a  one  dimensional  model. 
Substantial  computational  savings  are  achieved  when  compact  models  instead  of  high-fidelity 
models  are  used  for  synthetic  jets.  In  the  next  section  we  will  demonstrate  the  synthetic  jet 
arrays  simulations  for  two  practical  applications:  aerodynamic  control  of  airfoils  and  active  spot 
cooling  of  electronic  packages. 

4.2.7  Evaluation  of  CFD-ACE+  Mixed-Dimensional  Capabilities  with  Synthetic  Jets 

Microfabricated  arrays  of  synthetic  jets  are  being  explored  for  applications  in  aerodynamic  flow 
control,  in  cooling  of  electronic  packages,  and  in  mixing  in  microchemical  reactors. 
Computational  simulation  of  coupled  unsteady  fluid  mechanics  and  electromechanical  actuation 
of  a  single  synthetic  jet  can  be  performed  with  available  CFD  tools.  Modeling  of  flow  physics  of 
large  arrays  of  synthetic  jets  is  computationally  very  challenging.  At  CFDRC  we  have  explored 
two  complementary  computational  techniques:  a  three-dimensional  high-fidelity  model  for 
detailed  simulation  of  a  single  jet  (validated  with  measurements  results  from  GT),  and  a  reduced 
"single-cell"  model  of  a  jet  (described  above)  used  to  calibrate  the  compact  model. 

On  the  basis  of  the  3-D  simulations,  a  reduced  "single-cell"  model  of  synthetic  jet  has  been 
developed  at  CFDRC  (see  the  previous  section).  The  reduced  model  allows  for  simulation  of 
large  arrays  of  synthetic  jets,  for  example,  for  active  thermal  control  of  microelectronic  packages 
with  a  2D  array  of  synthetic  jets.  The  single-cell  models  have  been  inserted  into  3-D  mesh  in 
CFD-ACE,  used  for  high-fidelity  simulation  of  the  entire  array  of  synthetic  jets,  located  above  a 
high-power  electronic  module.  An  example  of  such  a  simulation  is  presented  in  Figure  99.  The 
objective  of  this  effort  is  to  analyze  active  control  of  synthetic  jets  for  spot  cooling  of  electronic 
packages. 

A  similar  technology  utilizing  Reduced  Models  of  Synthetic  Jets  will  be  used  within  the  other 
DARPA-funded  HERETIC  Program  (Heat  Removal  by  Thermo-Integrated  Circuits),  where 
CFDRC  will  support  modeling  and  design  of  Arrays  of  Synthetic  Microjets  for  active  cooling  of: 

-  VCSEL  Arrays  (Linear  and  2D)  -  Georgia  Tech  +  Honeywell 

-  Microprocessor  Chips  -  Georgia  Tech  +  Intel. 
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Figure  99.  3-D  simulation  of  7x7  array  of  synthetic  jets,  using  reduced  models:  Pressure 

Distribution. 

Reduced  Models  of  Synthetic  Jets  for  Aerospace  Applications 

Reduced  single-cell  model  of  synthetic  jet  has  been  used  at  CFDRC  for  simulation  of  arrays  of 
these  devices  in  the  aircraft  wing.  They  may  be  used  for  active  control  of  wing  aerodynamics  and 
other  possible  applications,  as  illustrated  in  Figure  100.  The  reduced  models  have  been  inserted 
into  3-D  high-resolution  mesh  in  CFD-ACE,  used  for  simulation  of  the  unsteady  aircraft  wing 
aerodynamics.  Figure  101  shows  the  CFD-ACE  simulation  results  of  airfoil  aerodynamics 
control  by  an  array  of  synthetic  jets. 

The  reduced  models  can  be  also  inserted  into  a  full  high-fidelity  simulation  of  an  airfoil,  to 
analyze  active  flow  control  (Figure  101).  A  large  array  of  synthetic  jets  is  positioned  on  an  airfoil 
and  demonstrated  on  active  aerodynamic  control  of  lift,  drag,  stall,  and  other  flow  characteristics. 
This  type  of  simulation  enables  "virtual  flight  control"  which  may  be  used  in  the  process  of 
design  of  Micro  Air  Vehicle  (Figure  102). 
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3-D  High-Fidelity  Models 


Compact/Reduced  Models 


System-Level  Simulation 
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Aerospace  Applications. 


Figure  101.  Results  of  CFD-ACE  simulations  at  CFDRC,  showing  the  use  of  reduced  models  of 
synthetic  jets  in  high-fidelity  3-D  simulations  of  airfoil:  overall  view,  leading  edge  and  jets,  syn-jets 

array  close-up. 


Active  jets  on  the  bottom  -  Take-Off  Active  jets  on  the  top  -  Landing 


Figure  102.  Micro  Air  Vehicle  (MAV),  and  its  virtual  flight  control,  vpith  3-D  CFD  simulations  and 

reduced  models  of  synthetic  jets. 

The  simulations  as  presented  above  demonstrate  CFD-ACE+  Mixed-Dimensionality  in  the 
design  with  applications  of  microdevices.  In  those  cases,  reduced  models  (one-dimensional,  or 
"single-cell"  models)  of  microdevices  were  embedded  in  the  high-fidelity  simulations  of  larger 
systems. 

4.2.8  Mixed-Dimensional  and  Mixed-Domain  Modeling  and  Experimental  Validation 

Micromachined  synthetic  jets  with  integrated  MEMS  modulators  have  been  fabricated  and 
characterized.  Several  methods  have  been  used  to  test  the  modulator  arrays:  visual  deflection 
test,  flow  modulation  test  by  smoke  visualization,  Schlieren  visualization,  pivot  tube,  x-wire 
anemometer,  and  particle  image  velocimetry  (PIV).  Most  of  these  methods  proved  to  be  very 
destructive  way  of  testing  the  devices  because  the  types  of  smoke  used  led  to  rapid  clogging  of 
the  jet  orifices. 

The  GT  team  devised  a  way  to  nondestructively  test  the  samples  by  viewing  the  cooling  effect  of 
our  modulated  synthetic  jets  upon  a  heated  plate  using  an  infrared  camera.  In  these  experiments, 
an  electric  heating  element  is  placed  on  the  copper  side  of  a  copper  clad  sheet  of  Kapton  (Figure 
103).  The  synthetic  jets  are  directed  at  the  heating  element,  and  an  infrared  camera  is  used  to 
view  the  temperature  distribution  of  the  heated  surface  cooled  by  the  synthetic  jets.  The 
temperature  distribution  on  the  heated  surface  varies  depending  upon  the  number  and  strength  of 
the  jets  impinging  upon  it. 

Figure  104  shows  sample  IR  camera  results  of  the  temperature  distribution  on  the  surface 
without  cooling  and  cooled  by  various  configurations  of  the  synthetic  jet  actuators.  The 
measured  data  are  compared  with  numerical  simulations  at  CFDRC,  including  use  of  reduced 
models  in  jet  arrays. 

Details  of  the  development  and  measurements  of  synthetic  jets  at  GT  were  described  in  Section 
4.2.2. 
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Figurel03.  Diagram  of  synthetic  jet  cooling  experiment. 
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Figure  104.  Sample  IR  camera  results  of  the  temperature  distribution  on  the  surface  without 
cooling  and  cooled  by  various  configurations  of  the  synthetic  jet  actuators. 


Mixed-Dimensional  and  Mixed-Domain  Simulation  of  Temperature  Control  by  Microjet 
Arrays 

At  CFDRC,  a  high-fidelity  3-D  model  has  been  built  to  run  a  set  of  simulations  corresponding  to 
the  GT  experiments,  as  described  in  the  previous  section.  The  CED-Micromesh  program  was 
used  to  build  a  full  3-D  model  of  the  jet  array  and  heated  substrate  (Figure  105). 


Figure  105.  A  high-fidelity  3-D  model  huilt  with  CFD-Micromesh,  to  run  a  set  of  simulations 
corresponding  to  the  Georgia  Tech  experiments  with  microjet  cooling. 


The  above  model  was  used  for  a  series  of  thermo-fluidic  simulations,  including  the  use  of 
reduced  models  developed  earlier  in  this  project.  The  full  model  consisted  also  of  an  array  of  2  x 
2  synthetic  jets,  a  cavity  as  a  jets  actuator,  and  a  heater  plate.  For  the  synthetic  jets,  the  single¬ 
cell  reduced  models  were  used,  within  the  thermal  enviroment  represented  by  a  full  3-D  model. 

A  fully  coupled  mixed-dimensional  and  mixed-domain  simulation  of  temperatutre  control  by 
microjet  arrays  was  performed.  As  a  result,  we  obtained  the  air  flow  pattern  caused  by  the 
synthetic  jets  (reduced  models)  and  temperature  changes  on  the  heater  plate  from  the  air-flow 
cooling  -  see  Figure  106.  The  results  were  also  compared  with  GT  measurements  for  models 
validation  and  calibration. 
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Figure  106,  Coupled  mixed-dimensional  and  mixed-domain  simulation  of  temperature  control  by 

microjet  arrays. 
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5  MODEL  REDUCTION  PROCEDURES 

5.1  Extraction  of  Lumped  Capacitances  from  3-D  Electrostatic  Simulations 

One  of  the  examples  of  coupled  problems  is  the  comb-drive  resonator,  used  in  MEMS 
accelerometers,  gyroscopes,  and  electronic  filters.  In  such  a  device,  electrostatic  actuation  is 
coupled  with  structural/dynamics  problems,  and  the  air  damping  is  very  important  too,  which  is  a 
fluidic-type  problem. 

As  one  of  the  steps  towards  automatic  generation  of  compact  models  of  comb  drives,  procedures 
have  been  developed  for  automatic  generation  of  compact  model  (extraction  of  lumped 
capacitances)  for  electrostatic  behavior  of  a  comb-drive  resonator,  using  results  of  high-fidelity 
3-D  electrostatic  simulations  with  FastBEM  from  CFDRC.  The  procedure  is  similar  to  the  one 
described  by  the  team  of  Analogy  Inc.  and  Robert  Bosch  GmbH  in  [Teegarden  1998]. 

The  procedure  is  illustrated  below. 
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♦  Micro-electro-mechanical  structure  design  imported  from  CAD/EDA 
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♦  3-D  model  of  the  structure 


♦  Extracted  capacitance  vs.  displacement 
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Rotor-Ground  Capacitance 


♦  Fitted  capacitance  functions  (Legendre  poiynomials) 

C  =  f(x,z) 

■O’ 

♦  Electrostatic  force  calculation 


Repeat  the  same  process  for  other  capacitances: 


The  fitted  capacitance  functions  may  be  used  directly  in  a  system-level  simulator,  like  Saber, 
where  they  can  be  coded  into  templates  of  the  devices,  using  the  MAST  language. 
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5.2  Generation  of  Equivalent  Circuit  Models  for  Microfluidics 


There  has  been  a  growing  interest  in  last  years  in  the  development  of  fluidic  microsystems 
containing  microfluidic  elements,  like  channels,  valves,  pumps,  etc.  Simulation  and  design  of 
microfluidic  systems  requires  various  level  models:  high-fidelity  (usually  three-dimensional,  or 
3-D)  models  [Yang  1998],  for  design  and  optimization  of  particular  elements  and  devices  (e.g. 
their  geometrical  structure),  as  well  as  system-level  models  allowing  for  VLSI-scale  simulation 
of  such  systems.  For  the  latter  purpose,  reduced  or  compact  models  are  necessary  to  make  such 
system  simulations  computationally  feasible  [Voigt  1998],  [Bourouina  1996]. 

In  this  project,  we  developed  a  design  methodology  and  practical  approach  for  generation  of 
compact  models  of  microfluidic  elements.  In  this  procedure  we  use  high-fidelity  3-D  simulations 
of  the  microfluidic  devices  to  extract  their  characteristics  for  compact  models,  and  subsequently, 
to  validate  the  compact  model  behavior  in  various  regimes  of  operation.  The  compact  models  are 
generated  automatically  in  the  formats  that  can  be  directly  used  in  SPICE  or  SABER. 

SPICE  and  SABER  Models  of  Tesla  Valve 

To  show  an  example  for  a  nonlinear  fluidic  device,  the  generation  of  compact  model  for  “Tesla 
valve”  is  described  in  detail  in  this  section.  Tesla  valve  is  one  of  the  no-moving-parts  (NMP) 
valves  used  in  micropumps  in  MEMS.  Its  principle  of  operation  is  based  on  the  rectification  of 
the  fluid  [Forster  1995].  For  the  same  pressure  drop,  the  flow  in  the  forward  direction  through 
the  valve  is  greater  than  the  flow  in  the  reverse  direction  (Figure  107).  Hence,  this  device  may  be 
considered  as  “fluidic  diode”. 


cross-section 


_ _ - — - 798  all  dimensions  in  microns 

Figure  107.  Geometry  and  dimensions  of  the  analyzed  Tesla  valve. 


In  order  to  generate  an  equivalent  circuit  model  (SPICE  or  SABER  format)  of  our  example 
microfluidic  device  -  Tesla  valve  -  we  need  to: 

1.  Build  a  3-D  model  and  computational  mesh  of  the  Tesla  valve. 

2.  Perform  a  full  3-D  simulation  of  fluid  flow  in  the  device,  both  steady  and  transient. 

3.  Perform  a  set  of  parametric  runs  to  extract  lumped  parameters  of  the  Flow-Pressure 
characteristics,  in  order  to  extract  parameters  for  a  compact  model  (equivalent  RL  circuit). 
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4.  Write  the  compact  model  in  the  format  compatible  with  SPICE  or  SABER. 


Performing  the  steps  1  to  3  involves  using  CFDRC's  new  software  CED-Micromesh  for 
automatic  building  of  3-D  solid  model  and  computational  mesh  (Section  3.1),  and  parametric 
computations  of  characteristics  with  CFD-ACE+  which  is  also  a  new  feature  added  in  this 
project  (Section  2.2).  An  example  unstructured  mesh  from  CFD-Micromesh  and  sample  results 
are  shown  in  Figure  108. 

(a) 


(c) 


Figure  108.  High-fidelity  simulations:  a)  3-D  unstructured  mesh,  h)  flow  velocity  (arrows)  and 
pressure  (colors)  inside  the  device  for  forward  flow,  c)  U- velocity  component  (along  the  X  axis). 

As  a  result  of  those  parametric  runs  of  3-D  simulations  in  CFD-ACE+,  we  obtain  the  static 
characteristics  of  flow  versus  pressure  drop,  as  shown  in  the  Figure  109  and  Figure  110  below. 
Open  circles  represent  data  obtained  from  3-D  simulations,  dark  dots  are  the  values  measured 
experimentally  at  the  University  of  Washington  (Prof.  Fred  Forster),  and  the  continuous  line  is 
an  analytical  fit  (obtained  in  Excel).  The  entire  static  characteristic  of  the  flow  in  both  directions 
is  presented  in  Figure  111. 
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Figure  109.  Forward  flow  static  characteristic  of  the  analyzed  Tesla  valve. 
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Figure  110.  Reverse  flow  static  characteristic  of  the  analyzed  Tesla  valve. 
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The  new  CFD-ACE+  version,  expanded  within  this  project,  can  generate  output  of  the 
parametric  runs  in  the  form  of  a  text  table,  as  in  the  example  below: 


***************iif****************************i*:************ 

Example  of  output  file:  Tesla_CHANNEL_OUTLET.dat 

★  ★★it****************-************************************* 


Mass  Flow 
P  mean 


Case 

1 

Iter. 

# 

100 

1.556542E+03 

6.869466E+03 

Case 

2 

Iter . 

# 

100 

2 .385028E+03 

1.288921E+04 

Case 

3 

Iter. 

# 

100 

3 .044021E+03 

1.857335E+04 

Case 

4 

Iter. 

# 

100 

3 .610177E+03 

2 .405173E+04 

Case 

5 

Iter . 

# 

100 

4.114709E+03 

2.938847E+04 

Case 

6 

Iter. 

# 

100 

4.574686E+03 

3 .461645E+04 

These  text  tables  can  be  easily  copied  and  pasted  into  any  convenient  post-processing  software. 
In  our  case,  we  used  Microsoft  Excel  to  draw  the  curves  and  fit  the  analytical  curves  (continuous 
lines  and  the  equations  in  Figure  109  and  Figure  110.  The  analytical  equations  are  needed  for  a 
SPICE  model. 

SPICE  model 

The  equivalent  circuit  model  of  the  Tesla  valve,  to  be  used  in  any  version  of  the  popular  SPICE 
simulator,  can  be  represented  in  the  form  of  a  series  connection  of  nonlinear  resistance  (R) 
element  and  linear  inductance  (L)  element,  as  illustrated  in  Figure  1 12. 
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Tesla  valve  model 


Figure  112.  Equivalent  RL  circuit  model  of  Tesla  valve. 

Most  SPICE  versions  do  not  accept  an  arbitrary  nonlinear  R  element,  so  it  should  be  replaced  by 
a  voltage-controlled  current  source  (VCCS),  using  the  analytical  relations  derived  in  Figure  109 
and  Figure  1 10. 


To  obtain  the  L  value  from  the  high-fidelity  simulations  of  any  microfluidic  device,  we  can  use 
the  relation  of  dynamic  behavior  of  fluid  pressure  vs.  flow  rate  and  its  electrical  equivalent, 
voltage  and  current  in  an  inductance: 


P  =  L^ 
dt 


V  = 


(5.1) 

The  L  value  can  be  easily  found  from  the  results  of  two  first  steps  of  transient  simulation,  as 
illustrated  in  the  figure  below: 


time  [s] 

Figure  113.  Transient  response  of  fluid  flow  to  step  pressure  change,  used  for  extraction  of  L  value 

in  the  equivalent  circuit  model. 
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The  final  form  of  the  SPICE  model  is  presented  in  the  listing  below.  The  nonlinear  resistance 
(coming  from  the  fitted  equations  in  Figure  109  and  Figure  110)  is  represented  by  voltage- 
controlled  current  source,  Gl,  and  the  linear  inductance  (derived  from  the  results  in  Figure  113) 
is  represented  by  LI. 


*  Equivalent-circuit  RL  model  of  fluid  flow  in  Tesla-45  valve 

*  inductance  -  fluid  inertia 

LI  1  2  1.413e-10 

*  VCCS  -  nonlinear  resistance 

Gl  2  0  Value  =  {0 . 5* (SGN(V(2 , 0) ) +1) *9222 . 7*PWR{V(2 , 0) , 0 . 6674)  + 

+  0.5* (SGN(V(2, 0) ) -1) *5960.1*PWR{V(2, 0) ,0.6042) ) 

V_dc  1  0  DC  OV 

***  STATIC  CHARACTERISTIC 

.DC  V_dc  -0.8  0.8  0.02 

.PRINT  DC  V(l)  I (Gl) 

.probe 

.END 


A  static  characteristic  of  the  Tesla  valve,  obtained  using  PSpice  with  the  input  deck  as  above,  is 
presented  in  Figure  114  below.  It  represents  very  well  the  measured  and  3-D-simulated 
characteristics  shown  earlier  in  Figure  111. 


-SOOmU  -MOOmU  0mU  4G0mU  SOQinU 

□  I(L1) 

U_dc 

Figure  114.  Static  characteristic  of  Tesla  valve,  calculated  with  SPICE. 

In  order  to  verify  dynamic  behavior  of  the  compact  model  in  SPICE,  a  transient  response  of  the 
Tesla  valve  flow  for  an  oscillatory  pressure  change  at  the  inlet  has  been  calculated  both  by  the  3- 
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D  simulation  using  CFD-ACE+,  and  using  the  equivalent  R-L  model  in  SPICE.  Example  results 
obtained  with  full  3-D  Navier-Stokes  solution  and  equivalent-circuit  RL  model  (solved  by 
PSpice),  for  periodic  pressure  oscillations  of  frequency  3085  Hz  (fluidic  resonance  of  Tesla 
valve  [Forster  1995])  and  amplitude  0.5  atm,  are  presented  in  Figure  115. 


3-D  CFD  Simulation  Results 


Equivalent  RL  Circuit  (SPICE) 


Tesla  45 
-  Transient 


Vol.  Flowluiymln] 
Pressure  *0.1  [Pa] 


Figure  115.  Transient  response  of  the  Tesla  valve  flow  for  an  oscillatory  pressure  change  at 
the  inlet,  calculated  by  the  3-D  simulation  with  CFD-ACE+,  and  using  the  equivalent  R-L 

model  in  SPICE. 


SABER  model 

A  behavioral  model  of  the  nonlinear  fluid  flow  in  Tesla  valve  has  been  also  implemented  for  the 
Saber  simulator  with  models  written  in  MAST,  an  analog  hardware  description  language  of 
Saber.  The  Saber-MAST  version  of  the  compact  flow  model  uses  directly  fluidic  quantities: 
pressure  drop  as  input  (across  variable)  and  fluid  flow  as  output  (through  variable). 

A  compact  model  of  the  Tesla  valve  for  Saber  simulator  can  be  generated  even  more  directly 
than  for  SPICE.  In  Saber,  a  predefined  template  for  fluidic  "Piecewise-Linear  Device" 
("q_p_pwl"  template)  can  be  used  for  the  non-linear  viscous  friction  model  (R).  The  tabularized 
pressure-flow  results  from  CED-ACE-i-  can  be  copied  directly  to  the  "q_p_pwl"  template.  A  test 
Saber  circuit  with  such  a  model  of  Tesla  valve,  and  also  a  hydraulic  pressure  source  available  in 
Saber,  is  shown  in  Figure  116.  The  listing  of  the  template  is  shown  under  the  figure.  The  table  of 
pressure-flow  data  in  the  template  is  generated  directly  in  that  format  by  CFD-ACE+  program,  as 
a  result  of  3-D  simulations. 
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Figure  116.  Test  circuit  in  Saber,  including  PWL  Tesla  valve  model  and  hydraulic  pressure  source 

for  testing. 

######################################################################## 

#  Saber  netlist  for  design  Tesla_static  # 

#  Created  by  the  Saber  Integration  Toolkit  4. 3 -2. 9  of  Analogy,  Inc.  # 
######################################################################## 

g  p  pwl .  Tesla_.45  pi :  Inlet  p2:  Outlet  =  data=  [ 

(-8.00E-01,  -5.12E+03), 

{-7.00E-01,  -4.75E+03), 

(-6.00E-01,  -4.35E+03), 

(-5.00E-01,  -3,92E+03) , 

(-4.00E-01,  -3.45E+03) , 

{-3.00E-01,  -2.92E+03) , 

{-2.00E-01,  -2.29E+03), 

(-9.97E-02,  -1.50E+03), 

(-5.98E-02,  -1.09E-f-03)  , 

(-2 .99E-02,  -6.96E+02) , 

(0,  0), 

(2.26E-02,  7.17E+02) , 

(4,30E-02,  1.13E+03), 

(6.87E-02,  1.56E+03) , 

(1.29E-01,  2 .39E4-03)  , 

{1.86E-01,  3 .04E+03) , 

(2.41E-01,  3.61E+03), 

(2.94E-01,  4.11E+03) , 

{3.46E-01,  4.57E+03) , 

{3.98E-01,  5,00E+03) , 

{4.48E-01,  5.40E+03) , 

(4.98E-01,  5.77E+03), 

(5.48E-01,  6.13E+03), 

(6.46E-01,  6.80E+03), 

{7.43E-01,  7.41E+03)] 

pressure. P_Drop  pi : Inlet  p2: Outlet  =  dc=0.0 
######################################################################## 


To  include  the  transient  effects  in  the  model  of  a  microchannel,  that  is  inertial  effects 
(inductance),  and  capacitance  due  to  the  fluid  compressibility,  we  can  use  the  "line_r"  template 
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available  in  Saber.  The  line_r  template  models  a  rigid  tube  used  for  the  transmission  of  hydraulic 
fluid. 

The  hydraulic  inductance  is  modeled  using  the  following  equation  [MAST,  1999]: 

AP  =  (rho  len/area)  •  dQ/dt 

and  the  hydraulic  capacitance  is  modeled  using  the  equation: 

Q  =  -(area-len/bulk)  •  dP/dt 

where  P  is  the  pressure,  Q  is  flow,  rho  is  mass  density  of  the  fluid,  len  is  the  channel  length,  area 
is  the  cross-section  area  of  the  channel,  and  bulk  is  the  fluid  bulk  modulus  value.  Figure  116 
shows  the  full  Saber  model  of  the  Tesla  valve,  including  the  viscous  resistance  part  (Tesla_R) 
and  inductance  (Tesla_L).  The  capacitance  due  to  the  fluid  compressibility  can  be  neglected  in 
most  cases  of  microfluidic  devices. 

A  static  characteristic  of  the  Tesla  valve,  obtained  with  Saber,  using  the  template  as  above,  is 
presented  in  Figure  117.  It  represents  very  well  the  measured  and  3-D-simulated  characteristics 
shown  earlier. 


Figure  117.  Static  characteristic  of  Tesla  valve,  obtained  in  SABER. 

The  compact  model  of  an  arbitrary  microfluidic  channel  presented  here,  in  the  form  of  nonlinear 
RL  circuit,  is  very  useful  to  a  wide  community  of  SPICE  users  who  are  very  often  also  designers 
of  microfluidic  MEMS  and  microsystems  in  the  domain  of  bio-medical  and  chemical 
engineering,  integrated  very  often  with  control  and  drive  electronics.  On  the  other  hand,  the 
behavioral  model  implemented  in  SABER  offers  more  direct  simulation  of  fluidic  phenomena, 
readily  available  to  designers  in  hydraulic  units,  without  additional  translation  to  electrical 
equivalents. 
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The  use  of  high-fidelity  3-D  numerical  simulations,  together  with  the  procedures  described 
above,  gives  a  possibility  to  obtain  very  useful  compact/behavioral  models  of  nonlinear 
microchannel  behavior  in  fluidic  microsystems,  for  wide  range  of  driving  conditions. 

Circuit  Modeling  of  Tesla-Micropump 

In  order  to  verify  to  performance  of  the  above  derived  equivalent  circuit  models  of  Tesla  valve, 
we  used  these  models  in  a  simulation  of  a  bigger  fluidic  system,  containing  8  Tesla  valves  and  a 
pumping  membrane.  A  3-D  model  of  the  Tesla-Micropump  with  eight  valves  and  a  snapshot  of 
flow  results  from  transient  3-D  CFD  simulation  is  shown  in  Figure  118  below.  The  flow  results 
show  a  time  instant  when  the  pumping  membrane  is  moving  up,  pulling  the  fluid  in.  It  can  be 
clearly  seen  how  the  left  branch  of  4  valves  is  “forward  biased”  and  the  right  branch  is  “reverse 
biased”. 


-5.58702 


Figure  118.  A  3-D  model  of  the  Tesla-Micropump  with  eight  valves  and  flow  results  from  transient 

3-D  CFD  simulation. 

Using  the  equivalent  circuit  model  of  a  single  Tesla  valve,  derived  in  the  previous  section,  we 
have  built  a  circuit  model  of  the  entire  pump  shown  in  the  figure  above.  First,  one  branch  of  the 
four  valves  was  modeled  by  means  of  the  four  elementary  circuits  connected  in  series,  and  the 
pressure  changes  forced  by  the  moving  membrane  were  represented  by  a  pressure  (voltage) 
source  -  see  Figure  119  below. 
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One  Tesla  valve 


Figure  119.  Equivalent  circuit  model  of  one  micropump  branch  with  four  Tesla  valves. 

The  above  model  was  written  as  an  input  deck  in  SPICE  format,  shown  in  the  listing  below.  The 
single  valve  model  is  coded  as  a  subcircuit  (SUBCKT)  called  Tesla  with  two  terminals  IN  and 
OUT,  and  then  it  is  used  many  times  in  the  SPICE  model  of  the  entire  pump. 


SPICE 

*  RL  model  of  fluid  flow  (uL/s)  in  4  Tesla  valves  in  series 

.SUBCKT  Tesla  IN  OUT 

*  inductance  -  fluid  inertia 
LI  IN  N1  5.087e-2 

*  VCCS  -  nonlinear  resistance  (calculated  from  uL/s) 

G1  N1  OUT  Valuer  {0 . 5* (SGN(V(N1 , OUT) ) +1) *0 . 0707*PWR {V(N1, OUT)  , 0 . 6674)  + 
+  0.5* (SGN{V(N1, OUT) ) -1 ) *0 . 0946 *PWR {V(Nl,OUT) ,0.6042) } 

.ENDS 

XI  1  2  Tesla 

X2  2  3  Tesla 

X3  3  4  Tesla 

X4  4  0  Tesla 

V_dc  1  0  DC  OV 

***  STATIC  CHARACTERISTIC 
.DC  V_dc  -8e4  8e4  2e3 


.END 


Comparison  of  transient  flow  results  obtained  both  from  high-fidelity  3-D  CFD  simulations  and 
with  the  equivalent  circuit  model  of  the  full  pump  (Figure  120)  are  shown  in  Figure  121. 
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Voltage 
=  Pressure 


Current  Source 
=  Forced  Flow 
by  Membrane 


Figure  120.  The  high-fidelity  CFD  results  (pressure  distribution)  in  the  micropump  with  eight 
Tesla  valves  and  its  full  equivalent  circuit  model. 


Results  of  3-D  CFD-ACE+  Simulations: 


time  [s] 


Equivalent  Circuit  (SPICE)  Results: 
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Figure  121,  Transient  flow  results  of  Tesla  micropump  obtained  from  high-fidelity  3-D  CFD 
simulations  and  with  the  equivalent  circuit  model  solved  by  SPICE. 
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The  validation  results  presented  above  show  the  compact  model  of  an  arbitrary  microfluidic 
channel  developed  in  this  project  and  presented  here  can  be  very  useful  to  a  wide  community  of 
designers  of  fluidic  microsystems  who  are  very  often  also  SPICE  users. 

5.3  Matrix  Reduction  Procedures  for  Electro-Mechanical  Microdevices 


5.3.1  UF  Model  Reduction  Techniques 

Within  this  project,  the  subcontractors  at  University  of  Florida  (UF)  worked  on  reduced  order 
models  for  electro-mechanical  devices,  built  on  the  basis  of  Finite  Element  Method  (FEM) 
models.  The  UF  team  implemented  the  model-reduction  code  using  both  the  Lanczos  algorithm 
(see  Lanczos  [1950],  Gallivan  &  Grimme  [1996]  and  Bahram  Nour-Omid  [1991])  and 
Wilson-Yuan-Dickens  (WYD)  algorithm  (see  Wilson,  Yuan  &  Dickens  [1982])  in  the 
CFD-ACEUN  code  of  CFDRC.  The  codes  were  tested  using  a  3-D  MEMS  beam  microstructure 
from  MIT.  We  also  studied  Lanczos  and  WYD  methods  for  higher  efficiency  and  accuracy. 

The  comparison  of  the  corresponding  first  several  modes  (up  to  seventh  mode)  from  the 
full-order  FEM  model  and  reduced-order  model  created  by  the  Lanczos  algorithm  was  very 
successful.  Also  the  step  responses  of  the  displacement  at  the  central  node  from  the  full-order 
model  and  the  reduced-order  model  (dofs  =  9)  agree  very  well.  The  transient  analysis  with  the 
reduced-order  model  takes  much  less  CPU  time  (0.033  sec.)  than  that  with  the  full-order  model 
(21.0  sec),  which  is  more  than  six  hundred  times  faster. 

Similarly  for  the  Wilson-Yuan-Dickens  (WYD)  algorithm,  the  comparison  of  the 
corresponding  modes  from  the  full-order  model  and  reduced-order  model  showed  that  the  first 
several  modes  (until  fourth)  found  out  by  the  WYD  reduction  method  agree  well  with  the  full 
order  modes.  Tests  with  serpentine  flexure  MEMS  structure  showed  that  transient  analysis  with 
reduced  model  can  be  up  to  10,000  times  faster  than  the  full  3-D  FEM  computation.  All  the 
details  and  results  are  presented  below. 

The  work  at  University  of  Florida  was  concentrated  on  developing  reduced-order  models  for 
MEMS  systems.  In  this  project,  we  have  implemented  our  model-reduction  code  using  both  the 
Lanczos  algorithm  (see  Lanczos  [1950],  Gallivan  &  Grimme  [1996]  and  BahramNour-Omid 
[1991])  and  Wilson-Yuan-Dickens  (WYD)  algorithm  (see  Wilson,  Yuan&  Dickens  [1982])  in 
the  CFD-ACEUN  code  of  CFDRC,  and  tested  the  codes  using  the  3-D  MEMS  beam 
microstructure  examples  in  our  previous  reports.  We  also  studied  Lanczos  and  WYD  methods 
for  higher  efficiency  and  accuracy.  In  Section  5.3. 1.1,  we  give  a  brief  description  of  our  work  in 
this  project .  In  Section  5.3. 1.2,  we  describe  the  validation  of  the  implementation  of  our  codes  in 
CFD-ACEUN.  In  Section  5.3.2,  we  present  the  detailed  implementation  of  our  model  reduction 
codes.  In  Section  5.3.3,  we  explain  our  research  on  Lanczos  and  WYD  algorithms  for  higher 
effciency. 

5.3.1.1  Research  Summary 

•  Validation  of  the  implemented  model  reduction  codes  by  testing  MEMS  models. 

•  Implementation  of  the  model  reduction  codes  into  CFD-ACEUN. 
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•  Development  of  the  corresponding  model  reduction  codes. 

•  Investigation  on  model  reduction  methods  (Lanczos  and  WYD)  for  MEMS  FEM  models. 

5.3.1.2  The  Validation  of  the  Implementation 

The  validation  of  the  implementation  of  Lanczos  algorithm 

In  this  subsection,  we  validate  our  implementation  of  Lanczos  algorithm  in  CFD-ACEUN  code 
by  testing  a  new  MEMS  beam  microstructure. 

Because  the  separated  mass  matrix  sky_  m  of  the  model  is  still  not  available  in  the  current  CFD- 
ACEUN  code,  we  have  to  read  in  the  model  information  from  input  data  files  directly,  which 
includes  the  stiffness  matrix  sky  _k,  the  mass  matrix  sky  _,  and  the  pointer  vector  ma  of  the 
original  model.  In  the  model  reduction  code,  the  reduced-order  model  will  be  calculated  from  the 
full-order  model  and  ready  for  the  transient  analysis  in  CFD-ACEUN  code.  In  the  updated  CFD- 
ACEUN  code  with  separated  mass  matrix  sky_m,  we  will  pass  the  three  one-dimensional  arrays 
sky_k,  sky_m,  and  ma  into  our  model  reduction  routines  as  arguments. 

Here  we  will  use  a  new  3-D  MEMS  beam  microstructure  (Figure  122)  to  validate  the 
implementation  of  our  model  reduction  code  (Lanczos  algorithm)  in  the  CFD-ACEUN  code.  We 
use  the  geometric  properties  and  the  material  properties  in  Hung  &  Senturia  [1999].  The 
geometric  properties  are  as  follows: 

Ig  =  200pm,  1^  =  50pm 

h  =  2.2pm,  w  =  40.0pm 

The  material  properties  are  as  follows: 

E  =  149GPa,  p  =  2330Kg/m^ 

The  transverse  step  loading,  at  central  point,  along  the  axis  is 

F=1000nN,  fort>0. 

We  use  twenty-eight  3-D  beam  elements  and  the  total  active  degrees  of  freedom  (dofs)  is  126. 
We  have  the  reduced-order  model  with  the  number  of  retained  degrees  of  freedom  nine.  The 
starting  vector  of  computing  the  reduced-order  model  is  the  random  vector.  The  comparison  of 
computed  eigenvalues  obtained  with  different  models  (full-order  and  reduced-order  model)  is 
shown  in  Table  5.1.  We  found  that  in  this  case  the  first  seven  eigenvalues  of  reduced-order 
model  (dofs  =  9)are  obtained  accurately.  The  Table  5.2  shows  the  comparison  of  CPU  time  to 
generate  different  reduced-order  models.  Speed  up  factors  are  calculated  with  respect  to  the 
reduced-order  model  (dofs  =  36).  We  found  that  the  time  taken  on  generating  the  reduced-order 
model  is  nearly  proportional  to  the  size  of  the  reduced-order  model. 
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The  Table  5.2  shows  the  comparison  of  CPU  time  to  generate  different  reduced-order  models. 
Speed  up  factors  are  calculated  with  respect  to  the  reduced-order  model  (dofs  =  36).  We  found 
that  the  time  taken  on  generating  the  reduced-order  model  is  nearly  proportional  to  the  size  of  the 
reduced-order  model. 
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Figure  122.  Finite-element  model  of  3-D  MEMS  beam  microstructure. 
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Table  5.1  Step  Response  of  the  3-D  elastic  beam  microstructure.  Comparison  of 
eigenvalues  from  full-order  model  to  reduced-order  model. 


Mode 

Full-order  (dofs  =  126) 

Reduced-order  (dofs  =  p) 

relative  error 

1 

■dhuohi 

1.7804  X  10“ 

0 

2 

3.5952  X  10“ 

0 

3 

5.8118  X  10“ 

0 

4 

wmEsasmm 

1.3905  X  10‘^ 

0.007% 

5 

1.4518  X  10‘" 

0.02% 

6 

2.5892  X  10'^ 

0.32% 

7 

0.56% 

8 

4.7345  X  10'^ 

42% 

9 

7.9633  X  10'" 

1.7243  X  10‘" 

53.8% 

Table  5.2  Step  response  of  the  3-D  elastic  beam  microstructure.  Comparison  of  CPU  time 
for  generating  the  reduced-order  model  with  different  size. 


Reduced-order  (dofs) 

run  time  (msec) 

speed  up  factor 

9 

25 

3.64 

18 

46 

1.98 

36 

91 

1 

The  comparison  of  the  corresponding  modes  from  the  full-order  mode!  and  reduced-order  model 
are  shown  in  Figure  123  to  Figure  132.  We  found  that  the  first  several  modes  (until  seventh)  are 
found  out  by  the  Lanczos  reduction  method  successfully,  which  verified  the  calculation  of 
eigenvalues  in  Table  5.1. 

The  comparison  of  the  step  response  of  displacement  at  the  central  node  along  the  axis  ^  from 
the  full-order  model  and  the  reduced-order  model  (dofs  =  9)  is  shown  in  Figure  133.  The  Figure 
134  is  the  frequency  components  of  the  time-domain  response  by  Fast  Fourier  Transform  (FFT), 
which  shows  the  comparison  of  the  corresponding  frequency  response  from  the  full-order  model 
and  from  the  reduced-order  model  (dofs  =  9).  The  corresponding  velocity  comparisons  are  given 
in  Figure  135  and  Figure  136.  The  results  from  the  full-order  model  and  the  reduced-order  model 
(dofs  =  9)  agree  very  well.  On  the  other  hand,  the  transient  analysis  with  the  reduced-order 
model  takes  much  less  CPU  time  (0.033  sec.)  than  that  with  the  full-order  model  (21.0  sec), 
which  is  more  than  six  hundred  times  faster. 

From  the  above  numerical  results,  we  show  that  our  implementation  of  the  Lanczos  algorithm  in 
the  CFD-ACEUN  code  was  successful.  We  also  establish  the  efficiency  of  the  reduced-order 
model  by  Lanczos  algorithm. 
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Figure  123.  mode  in  full-order  model. 


Figure  124.  1®‘  mode  in  reduced-order  model 
(dofs  =9). 
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Figure  125.  2”**  mode  in  full-order  model.  Figure  126.  2"'*  mode  in  reduced-order  model 

(dofs  =9). 


Figure  127.  3"**  mode  in  full-order  model.  Figure  128.  3"*  mode  in  reduced-order  model 

(dofs  =9). 


Figure  129.  4***  mode  in  full-order  model.  Figure  130.  4***  mode  in  reduced-order  model 

(dofs  =9). 
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(dofs  =9). 


The  Validation  of  the  implementation  of  WYD  algorithm 

In  this  subsection,  we  will  validate  our  implementation  of  WYD  algorithm  in  CFD-ACEUN  code 
by  testing  the  MEMS  beam  microstructure  (Figure  122)  described  in  the  previous  subsection. 

We  use  twenty-eight  3-D  beam  elements  and  the  total  active  degrees  of  freedom  (dofs)  is  126. 
We  have  the  reduced-order  model  with  the  number  of  retained  degrees  of  freedom  nine.  The 
starting  vector  of  computing  the  reduced-order  model  is  the  static  solution. 

The  comparison  of  computed  eigenvalues  obtained  with  different  models  (full-order  and 
reduced-order  model)  is  shown  in  Table  5.3.  We  found  that  in  this  case  for  the  simple  constant 
concentrated  force  WYD  method  does  not  skip  any  of  the  first  fourth  modes.  The  mode 
participation  factor  shows  that  the  first  several  modes  are  all  important.  No  mode  is  orthogonal  to 
the  applied  force  as  in  our  previous  crab-leg  model  (see  Fedder  [1994]).  We  have  the  same  above 
observation  as  in  both  WYD  algorithm  and  Lanczos  algorithm  with  static  solution  as  starting 
vector,  however,  it  seems  WYD  algorithm  may  pick  up  only  one  mode  among  some  modes  with 
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close  eigenvalues.  The  comparison  of  the  corresponding  modes  from  the  full-order  model  and 
reduced-order  model  shows  that  the  first  several  modes  (until  fourth)  found  out  by  the  WYD 
reduction  method  agree  well  with  the  full  order  modes,  which  verified  the  calculation  of 
eigenvalues  in  Table  5.3.  The  comparison  of  the  step  response  of  displacement  at  the  central 
node  along  the  axis  ^3  from  the  full-order  model  and  the  reduced-order  model  (dofs  =  9)  is 
shown  in  Figure  137.  The  Figure  138  is  the  frequency  components  of  the  time-domain  response 
by  Fast  Fourier  Transform  (FFT),  which  shows  the  comparison  of  the  corresponding  frequency 
response  from  the  full-order  model  and  from  the  reduced-order  model  (dofs  =  9).  The  results 
from  the  full-order  model  and  the  reduced-order  model  (dofs  =  9)  match  very  well.  On  the  other 
hand,  the  transient  analysis  with  the  reduced-order  model  takes  much  less  CPU  time  (0.033  sec.) 
than  that  with  the  full-order  model  (23.0  sec),  which  is  more  than  six  hundred  times  faster. 


Figure  133.  Step  responses  of  displacement  at  central  point  from  full 
order  model  and  reduced-order  model  (dofs  =  9). 


From  the  above  numerical  results,  we  show  that  our  implementation  of  the  WYD  algorithm  in 
the  CFD-ACEUN  code  was  successful.  We  also  established  the  efficiency  of  the  reduced-order 
model  by  WYD  algorithm 

5.3.2  Implementation  of  Model  Reduction  Algorithms  in  CFD-ACEUN 
Implementation  of  the  Lanczos  model  reduction  code. 

Next,  we  will  describe  the  implementation  of  our  model  reduction  code  into  CFD-ACEUN.  We 
modified  the  make  file  of  CFD-ACEUN  to  accommodate  the  routines  related  to  our  model 


118 


reduction  code  which  includes  the  LANZ  package,  and  the  interface  between  LANZ  and  CFD- 
ACEUN.In  the  CFD-ACEUN,  the  routine  fldprs_fem.f  will  reads  in  the  finite  element  model 
from  the  preprocessor,  and  then  the  routine  startfem.f  set  up  the  element  properties  for  solving 
the  problem.  Then  the  element  consistency  is  checked  in  routine  check_elements:f.  After  that, 
the  routine  solved.f  will  assemble  the  element  stiffness  matrix  and  element  mass  matrix  and 
apply  the  corresponding  boundary  conditions,  and  form  the  global  matrix  and  residual  force.  The 
skyline  solver  will  be  used  to  solve  the  linear  algebra  system.  Finally,  the  routine 
cfdview_output:f  will  create  unstructured  CFD-VIEW  files  for  postprocessors.  We  will  add  the 
model  reduction  code  in  the  routine  solved.f. 

To  call  the  model  reduction  code,  we  are  adding  the  subprogram  model_reduction  inside  the 
routine  solved.f.  In  the  subprogram  modeLreduction,  the  routine  lanzjo  and  the  routine 
lanz_param  will  prepare  the  necessary  data  of  the  original  model  for  LANZ.  The  routine 
lanz_driver  will  then  call  LANZ  package  for  calculating  reduced  model.  All  the  above  three 
routines  will  be  put  in  one  module  file,  lanz_mod:f,  which  plays  the  role  of  the  interface  between 
LANZ  and  ACEUN.  Another  module  lanz_variable_mod:f  will  declare  all  variables  used  in  the 
reduced  model  code,  which  may  be  combined  with  the  ACEUN  variable  declaration  module  file 
later  on. 


Figure  134.  Frequency  responses  at  central  point  from  full-order  model 
and  reduced-order  model  (dofs  =  9). 


In  the  LANZ  software  package,  the  routine  lanz.fileS.f  is  the  main  subroutine  for  the  Lanczos 
algorithm,  which  was  coded  by  Fortran??.  The  routine  lanz.file4.c  is  a  C  routine,  which  allocates 
and  deallocates  memory  inside  LANZ  package.  We  have  modified  the  LANZ  in  order  to  use 


119 


Fortran90  function  ALLOCATE  to  allocate  the  arrays  used  in  LANZ.  Moreover,  we  have  coded 
the  routine  lanz.fileS.f  inside  the  LANZ  package,  to  obtain  the  reduced  model,  including  the 
reduced  matrices,  reduced  force  vector,  and  the  reduced  initial  conditions. 

In  our  reduced  model  eode,  the  stiffness  matrix  sky_k  and  the  mass  matrix  sky_m,  and  the 
pointer  vector  ma,  are  read  in  from  ACEUN.  The  matrices  sky_k  and  sky_m  are  stored  as  one¬ 
dimensional  arrays  in  ACEUN,  and  the  vector  ma  points  to  the  loeation  of  the  main  diagonal 
terms  in  the  skyline  storage.  Until  now  the  separated  mass  matrix  sky_m  ,  however,  is  not 
available  in  ACEUN.  We  appreciate  the  sky_m  can  be  given  in  the  updated  CFD-ACEUN  soon. 
For  the  convenience  of  testing  our  code,  we  also  wish  the  beam  element  will  be  added  in  the 
element  library  of  CFD-ACEUN.Once  the  stiffness  matrix  sky_k,  mass  matrix  sky_m,  the 
pointer  vector  ma,  and  the  initial  condition  of  the  original  model  are  delivered  from  ACEUN,  the 
model  reduction  part  in  ACEUN  will  calculate  the  Lanczos  vectors  and  the  reduced  model, 
including  a  symmetric  tridiagonal  matrix  with  the  required  reduced  size  and  the  corresponding 
reduced  initial  condition,  then  the  dynamic  response  of  the  reduced  model  can  be  computed. 

With  the  response  of  the  reduced  model  and  the  calculated  Lanczos  vectors,  the  dynamic 
response  of  the  original  model  will  then  be  recovered.  We  will  report  the  testing  of  the 
implementation  of  the  reduced  model  code  (Lanczos)  in  ACEUN. 


Figure  135.  Step  responses  of  velocity  at  central  point  from  fuli-order 
model  and  reduced-order  model  (dofs  =  9). 
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Implementation  of  the  WYD  model  reduction  code. 

Here,  we  will  further  explain  each  routine  in  our  model  reduction  codes.  Inside  the  CFD- 
ACEUN  routine  solved.!,  we  have  coded  a  subroutine  model_reduction  for  generating  the 
reduced  model  (WYD),  which  calls  subroutine  WYD_driver.  WYD_driver  is  put  in  the  module 
_le  WYD_mod:f.  In  addition  to  the  Gram-Schmidt  or-thogonalization  process.  Inside  subroutine 
WYD_driver,  we  have  implemented  selective  reorthogonalization  to  ensure  the  linear 
independence  among  the  generated  WYD  Ritz  vectors.  All  variables  used  in  the  model  reduction 
code  (WYD)  are  put  in  another  module  file,  wyd_variable_mod.f.  We  are  describing  the 
arguments  of  the  above  subprograms  and  module  files  as  follows: 

WYD_driver(reduce_active_dof,  ma;  sky_k, 

sky_m,  WYD_vector,  sky_rhs,  reduce_rhs, 
reduced_k,  reduced_m,  ini_disp, 
ini_vel,  reduce_ini_disp,  reduce_ini_vel) 
reduce_active_dof :  input,  total  equation  number  of  the  reduced  model, 
ma  :  input,  vector  of  the  stiffness  matrix, 
sky_k  :  input,  stiffness  matrix  of  the  original  structure, 
sky_m  ;  input,  mass  matrix  of  the  original  structure, 

WYD_vector :  output,  WYD  vectors  for  recovering  the  response, 
sky_rhs  :  input,  the  right  hand  side  of  the  original  model, 
reduce_rhs  :  output,  spatial  external  force  of  the  reduced  mode,l 
reduced_k  :  output,  stiffness  matrix  of  the  reduced  model, 
reduced_m  :  output,  mass  matrix  of  the  reduced  model, 
ini_disp  :  initial  displacement  vector  of  the  original  model, 
ini_vel :  input,  initial  velocity  vector  of  the  original  model, 
reduce_ini_disp  :  output,  initial  displacement  vector  of  the  reduced  model, 
reduce_ini_vel :  output,  initial  velocity  vector  of  the  reduced  model. 

matrix_multip_sky(sky_index,  mass,  x,  v) 

skyjndex  :  input,  point  of  vector  of  the  skyline  matrix  mass, 
mass  ;  input,  mass  matrix  to  perform  multiplication  wrt  Ritz  vectors, 

X  :  input,  Ritz  vector, 

V  :  output,  product  of  mass  matrix  and  a  Ritz  vector. 

matrix_multip(reduced_dof,  x,  a,  result_a) 

reduced_dof :  input,  total  equation  number  of  reduced  model, 

X  :  input,  matrix  to  perform  multiplication  wrt  Ritz  vectors, 
a  :  input,  Ritz  vector, 

V  :  output,  product  of  matrix  and  a  Ritz  vector 

matrix_product(reduced_dof,  skyjndex,  a,  x,  redu_a) 

reduced_dof :  input,  total  equation  number  of  reduced  model 
skyjndex  :  input,  point  vector  of  the  skyline  matrix  mass 
X  :  input,  matrix  to  perform  multiplication  wrt  Ritz  vectors, 
a  :  input,  Ritz  vector. 
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redu_a  ;  output,  product  of  a  Ritz  vector,  matrix  and  Ritz  vector, 

m_normal(sky_index,  mass,  force,  n_vector) 

skyjndex  :  input,  point  vector  of  the  skyline  matrix  mass, 
mass  :  input,  mass  matrix  a  Ritz  vector  will  be  normalized  to, 
force  :  input,  a  Ritz  vector  will  be  normalized, 
n_vector  :  output,  a  normalized  Ritz  vector. 


Frequecicjy  response  ot  full  order  model 


Figure  136.  Frequency  responses  of  velocity  at  central  point  from  full- 
order  model  and  reduced-order  model  (dofs  =  9). 
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Figure  137.  Step  responses  of  displacement  at  central  point  from  full- 
order  model  and  reduced-order  model  (dofs  =  9). 

5.3.3  Research  on  the  Model  Reduction  algorithms. 

In  this  section,  we  present  our  study  on  the  Lanczos  Algorithm  with  static  starting  vector  in 
Subsection  5.3.3. 1.  In  Subsection  5.3.3.2,  we  describe  our  study  on  the  Loss  of  orthogonality  and 
reorthogonalization  in  WYD  algorithm  (see  Leger  [1986]  and  Paige  [1972]).  In  Subsection 
5. 3. 3.3,  we  test  the  Karhunen-Loeve  decomposition  model  reduction  method  used  by  Hung  & 
Senturia  [1999].  Since  a  few  full  order  model  runs  are  needed  to  generate  reduced-order  model 
by  Karhunen-Loeve  decomposition,  the  overall  computational  cost  is  much  higher  than  our 
WYD  model  reduction  method  with  the  same  response  accuracy.  We  also  studied  error 
estimation  in  the  WYD  algorithm  based  on  the  loading  representation.  The  participation  factors 
given  in  Table  5.3  clearly  show  which  modes  are  important  and  which  modes  may  be  skipped  by 
the  WYD  algorithm. 
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Figure  138.  Frequency  responses  at  central  point  form  full-order  model  and  reduced  - 

order  model  (dofs  =  9). 

5.3.3.1  Lanczos  Algorithm  with  Static  Starting  Vector 

Unlike  WYD  method,  the  starting  vector  in  the  Lanczos  method  is  usually  chosen  at  random, 
ignoring  the  important  information  specific  to  the  dynamic  problem.  The  Lanczos  method  is 
expected  to  be  more  effective  if  the  spatial  distribution  of  the  dynamic  load  is  used  to  initiate  the 
recurrence  relationship  in  Lanczos  vector  generation  process.  The  numerical  results  show  the 
Lanczos  method  with  static  starting  vector  performs  much  better  than  random  starting  vector.  We 
have  coded  it  as  one  option  when  using  Lanczos  method  in  the  model  reduction  code.  Here  we 
use  the  same  example  of  the  crab-leg  flexure. 

We  still  use  3-D  Euler-Bemoulli  beam  elements  to  model  the  crab-leg  exure.  We  use  the  same 
geometric  properties  and  the  material  properties  from  Fedder  [1994]  as  well.  We  use  twenty  3-D 
beam  elements  and  the  total  active  degrees  of  freedom  (dofs)  is  78.  We  have  two  cases  of 
reduced-order  model,  i.e.,  the  number  of  retained  degrees  of  freedom  are  nine  and  five.  The 
comparison  of  computed  eigenvalues  obtained  with  different  models  (full-order  and  reduced- 
order  model)  is  shown  in  Table  5.4. 

If  we  use  the  random  starting  vector,  we  find  that  the  first  three  of  reduced-order  model  (dofs  = 
9)  and  two  eigenvalues  of  reduced-order  model  (dofs  =  5)  are  obtained  accurately.  Now  we  use 
the  static  solution  as  the  starting  vector,  we  find  that  the  eigenmodes  obtained  are  not  necessarily 
the  first  modes  in  full  model.  For  example,  in  the  case  of  reduced-order  model  (dofs=5),  the  first 
eigenmode  corresponds  to  the  first  eigenmode  in  the  original  model,  the  second  eigenmode 
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corresponds  to  the  sixth,  the  third  eigenmode  corresponds  to  the  seventh  eigenmode  (see  Table 
5.4). 

The  comparison  of  the  step  response  at  the  central  node  along  the  axis  _from  the  full-order 
model  and  the  reduced-order  model  (dofs  =  9)  is  shown  in  Figure  133.  Figure  134  is  the 
frequency  components  of  the  time-domain  response  by  Fast  Fourier  Transform  (FFT),  which 
shows  the  comparison  of  the  corresponding  frequency  response  from  the  full-order  model  and 
from  the  reduced-order  model  (dofs  =  9). 

From  the  above  figures,  we  find  that  much  better  results  are  obtained  with  the  static  starting 
vector.  It  is  observed  when  the  size  of  the  reduced-order  model  is  5,  the  responses  in  both  time 
domain  and  frequency  domain  already  agree  very  well,  comparing  to  the  full-order  model.  The 
Lanczos  method  with  static  solution  as  the  starting  vector  is  much  better  than  with  random 
starting  vector.  It  is  because  the  eigenmodes  picked  here  in  Lanczos  method  are  the  most 
important  modes  for  dynamic  response,  not  the  first  modes  of  the  original  model.  Therefore,  we 
strongly  recommend  to  use  Lanczos  method  with  the  static  starting  vector. 

On  the  other  hand,  Lanczos  method  can  achieve  the  same  good  results  as  WYD  method,  but  with 
numerically  a  more  efficient  reduced  model,  i.e.,  a  symmetric  tridiagonal  matrix  Tr  with  the 
reduced  size: 


LY  +  IY  =  L  (5.2) 

In  WYD  method,  we  obtained  a  full  symmetric  matrix  Kr  with  the  reduced  size  in  the  reduced 
model: 


IY  +  KrY  =  L 


where  /  is  the  identity  matrix  with  reduced  size.  Obviously,  Eq.  (5.2)  will  be  solved  more  quickly 
than  Eq.  (5.3). 
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Figure  139.  Relative  step  response  error  of  displacement  at  central  point  bewteen  full- 
order  model  and  reduced-order  model  (dofs  =  9). 

5.3.3.2  Loss  of  Orthogonality  and  Reorthogonalization  in  WYD  Algorithm 

The  selective  reorthogonalization  is  used  to  maintain  the  necessary  orthogonality  among  the  trial 
vectors  in  WYD  method.  In  this  subsection,  we  will  show  that  random  starting  vector  is  better 
than  the  static  solution  as  starting  vector  in  terms  of  orthogonality  among  the  trial  vectors.  Loss 
of  orthogonality  happens  until  the  spectral  content  of  the  starting  vector  is  exhausted  if  we  use 
random  starting  vector,  and  we  can  obtain  a  set  of  independent  vectors  which  has  the  same 
number  as  the  original  order  in  the  iterative  process  even  without  reorthogonalization.  Instead  of 
using  random  starting  vector,  if  we  use  static  solution  as  starting  vector,  the  loss  of  orthogonality 
will  happen  when  the  number  of  generated  WYD  vectors  is  greater  than  half  of  the  number  of 
the  original  order  without  reorthogonalization.  Even  with  reorthogonalization,  the  independency 
of  generated  WYD  vectors  is  sensitive  with  respect  to  the  reorthogonal  tolerance. 

We  still  use  3-D  Euler-Bemoulli  beam  elements  to  model  the  crab-leg  flexure.  We  use  the  same 
geometric  properties  and  the  material  properties  from  Fedder[1994]  as  well.  We  use  twenty  3-D 
beam  elements  and  the  total  active  degrees  of  freedom  (dofs)  is  78.  If  we  use  random  starting 
vector  the  vector  basis  can  be  generated  even  without  reorthogonalization  until  the  spectral 
content  of  the  starting  vector  is  exhausted.  The  reduced  number  of  dofs  can  keep  going  to  the 
number  of  72  without  reorthogonalization.  If  we  use  static  solution  as  starting  vector,  loss  of 
orthogonality  will  happen  when  the  number  of  generated  WYD  vectors  is  greater  than  60  without 
reorthogonalization.  After  reorthogonalization,  the  number  of  independent  WYD  vectors  can 
increase  to  74. 
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Figure  140  shows  the  global  orthogonality  obtained  for  different  starting  vector  without 
reorthogonalization  from  the  original  WYD  algorithm  in  our  code.  It  clearly  shows  that  random 
starting  vector  will  not  loss  orthogonality  among  trial  vectors  till  the  end  of  the  generation 
process. 

A  random  starting  vector  is  usually  chosen  hoping  that  all  required  eigenvectors  will  be 
represented.  Another  observation  is  loss  of  orthogonality  still  may  happen  even  after  a  few 
number  of  reorthogonalization  iterations.  That  means  we  do  not  give  a  small  enough  tolerance  as 
orthogonal  criteria.  This  can  be  solved  if  we  decrease  the  tolerance,  from  l.OE-16  to  l.OE-20,  for 
instance.  However,  it  will  cost  more  computational  effort.  Figure  141  shows  the  global 
orthogonality  obtained  for  static  solution  as  starting  vector  before  and  after  reorthogonalization 
from  the  original  WYD  algorithm  in  our  code.  It  shows  that  even  after  reorthogonalization  it  still 
may  lose  orthogonality  among  trial  vectors  in  the  end  of  generation  process. 

Our  results  show  that  random  starting  vector  is  good  to  generate  almost  the  full  set  of  WYD 
vector  basis  even  without  reorthogonalization.  Therefore  to  study  the  sensitivity  of  the 
orthogonal  degree  among  the  generated  WYD  vectors  with  respect  to  the  orthogonal  tolerance, 
we  use  static  solution  as  starting  vector. 

In  Figure  142,  we  compared  the  relationship  between  number  of  iterations  and  the  number  of 
WYD  vectors  for  different  tolerance. 

Our  result  shows  that  to  generate  a  set  of  well  independent  WYD  vectors,  we  need  to  be  very 
careful  in  choosing  the  tolerance.  If  the  tolerance  is  not  small  enough,  the  generated  vectors  will 
not  be  exactly  independent  (here  not  exactly  independent  means  not  independent  enough  to 
ensure  the  following  WYD  vectors  being  generated  to  be  orthogonal  with  respect  to  vectors 
already  been  generated).  Although  it  still  can  keep  a  satisfactory  degree  of  orthogonality  for 
generating  a  few  number  of  WYD  vectors,  it  will  cost  more  computational  effort  if  we  want  to 
generate  more  independent  WYD  vectors.  The  number  of  reorthogonal  iterations  will  have  a 
great  increase  in  order  to  make  the  vector  being  generated  orthogonal  with  respect  to  the 
previous  generated  vectors.  On  the  other  hand,  if  the  tolerance  is  too  small,  it  will  obviously 
increase  the  number  of  reorthogonal  iterations.  For  the  crab-leg  model  Figure  142  shows  that 
there  is  an  optimal  tolerance  value  TOL=1.0E-20  because  of  its  less  overall  computational  effort 
(less  number  of  reorthogonal  iterations).  In  our  WYD  model  reduction  code,  we  leave  this  as  an 
option  for  users. 
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Orttio^fonal  Vector  Basts  —  crab^leg  beam  model;  orlglnai  DOF  =  78 


Number  of  Calculated  WYD  Vectors 

Figure  140.  The  global  orthogonality  level  maintained  by  the  original  WYD  algorithm  for 

different  starting  vectors. 

crab_ieg  beam  model  (static  solution  as  starting  vector);  original  CX>F  =  7& 


Number  of  Calculated  WYD  Vectors 

Figure  141.  The  global  orthogonality  level  maintained  by  the  original  WYD  algorithm  for 

static  solution  as  starting  vector. 
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Sen^livtly  of  incfepl  WYD  vector  wrt  reortfwQ  tolerance;  original  DOF  =  78 


0  10  20  30  40  50  60  70  80 

Number  of  WYD  Vectors 

Figure  142.  Sensitivity  of  independent  WYD  vector  with  respect  to  reorthogonal  tolerance. 


Table  5.3.  Step  response  of  the  3-D  elastic  beam  microstructure.  Comparison  of 
eigenvalues  from  full-order  model  to  reduced-order  model. 
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7 

3.2409  X  10'" 

3.2409  X  10'^ 

4.8759  X  10'^ 

8 

4.7345  X  10'" 

4.7345  X  10‘^ 

2.1281  X  10'^ 

9 

7.9633  X  10*^ 

7.9635  X  10‘^ 

2.2695  X  10  " 

10 

8.7063  X  10'" 

8.7222  X  10*^ 

5.7315  X  10'^ 

11 

1.2039  X  10'" 

1.2041  X  10*^ 

6.4250  X  10’'' 

12 

1.4504  X  10'" 

1.5566  X  10'" 

6.1548  X  10"^ 

13 

1.8232  X  10'^ 

2.6336  X  10'-^ 

1.6704  X  10^' 

14 

2.5171  X  10'" 

3.8518  X  10'^ 

2.7047  X  10'^ 

15 

2.7105  X  10'^ 

1.2596  X  10‘^ 

2.3274  X  10'^ 

16 

3.2016  X  10‘^ 

6.3128  X  lO'^ 

5.9028  X  10'^ 

Table  5.4  Step  response  of  the  crab-leg  flexure.  Comparison  of  eigenvalues  form  full-order 


model  to  reduced-order  model. 


mode 

full-order  (dofs  =  78) 

WYD  (dofs  =  9) 

reduced-order  (dofs  =  5) 

1 

1.3869  X  lO"* 

1.3869  X  10^ 

1.3869  X  lO"' 

2 

4.6851  X  109 

4.6848  X  lO"' 

9.8316  X  10*^ 

3 

5.3343  X  10"* 

9.8316  X  lO"* 

1.7985  X  10'^ 

4 

8.1038  X  lO'^ 

1.7985  X  10**^ 

8.1502  X  10'" 

5 

9.1761  X  10"' 

8.1501  X  10"^ 

7.5029  X  10‘^ 

6 

9.8316  X  lO"" 

3.3673  X  10‘‘ 

7 

1.7985  X  10'" 

6.0128x10"  n 

8 

1.9378  X  10“^ 

1.0280  X  lO*'' 

9 

2.4627  X  10'^ 

1.4154  X  10‘^ 

5.3.3.3  Comparison  of  WYD  Algorithm  and  Kahrunen-Loeve  Decomposition 

Karhunen-Loeve  (KL)  decompositions  were  originally  used  to  generate  basis  functions  for 
turbulence  problems  in  fluid  mechanics,  fluid  structure  interactions  in  aerodynamical  systems, 
and  in  chemical  engineering  systems.  KL  procedure  uses  the  so-called  snapshot  method,  where 
the  problem  of  obtaining  eigenmodes  of  a  large  system  reduces  to  solving  eigenmodes  of 
matrices  of  order  of  only  100.  Second,  the  method  always  produces  real,  optimal  modes 
regardless  of  the  damping  characteristics  of  the  system  under  consideration.  Third,  the  method  is 
a  direct  response  approach  that  does  not  require  a  dynamic  model  describing  the  system.  In  this 
section,  we  present  some  test  results  using  the  MEMS  microstructure  by  KL  decomposition.  Our 
simulation  procedure  using  MATLAB  is  as  follows: 

•  Construct  the  state  space  form  of  the  full-order  dynamic  model. 

•  Run  a  few  time  steps  to  sample  the  state  variables  at  a  series  of  different  times. 

•  Do  single  value  decomposition  for  the  state  matrix  consists  of  selected  state  vectors. 

•  Pick  up  basis  vectors  as  Ritz  vectors  to  construct  reduced  order  model. 

The  comparison  of  the  relative  response  error  of  displacement  and  velocity  at  the  central  node 
along  the  axis  ^  from  the  WYD  algorithm  and  the  Karhunen-Loeve  decomposition  is  shown  in 
Figure  139  -  Figure  144.  Table  5.5  shows  the  comparison  of  computime  for  generating  the 
reduced  order  model  using  WYD  algorithm  and  Karhunen-Loeve  decomposition.  It  was  clearly 
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shown  that  with  the  same  order  of  response  accuracy,  Karhunen-Loeve  decomposition  has  a 
much  higher  computational  cost  than  WYD  algorithm. 


Table  5.5.  Step  response  of  the  3-D  elastic  beam  microstructure.  Comparion  of  CPU  time 
for  generating  the  reduced-order  model  with  different  algorithms. 


Reduced-order  (dofs) 

Run  time  of  WYD  (sec) 

Run  time  of  KL  (sec) 

9 

0.12 

22.6 

18 

0.30 

23.1 

36 

0.77 

22.8 

0  0.5  1  t.S  2  2.5  3  3.5  A 

*  X 10*’ 

Figure  143.  Relative  step  response  of  error  of  displacement  at  central  point  between  full- 
order  model  and  reduced-order  model  (dofs  =  9). 
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Figure  144.  Relative  step  response  error  of  velocity  at  central  point  between  full-order 

model  and  reduced-order  model  (dofs  =  9). 

5.4  Implementation  of  MIT  Model  Reduction  Procedure 

In  the  frame  of  this  project,  we  have  also  implemented  at  CFDRC  the  procedures  of  FEM  model 
reduction  based  on  Arnold!  method  for  electro-mechanical  microdevices,  developed  at  Prof. 
Jacob  White’s  group  at  MIT  within  the  Composite  CAD  Program.  The  codes  have  been 
implemented  and  installed  at  CFDRC  and  coupled  with  CFD-ACE+  code  by  Deepak 
Ramaswamy,  PhD  student  from  MIT,  during  his  summer  internship  at  CFDRC  in  Huntsville  in 
August-September  2000.  The  new  Amoldi-based  model  reduction  procedures  were  demonstrated 
and  tested  in  coupling  with  CFD-ACE-t-  code. 


The  Model  Reduction  (MR)  Procedure  steps  are: 

1.  Build  a  full  3-D  FEM  model  of  the  electro-mechanical  structure. 

2.  Perform  one  steady-state  coupled-domain  simulation  (structures  +  electrostatics)  with  the 
full  FEM  model,  for  the  applied  bias  Vo,  using  CFD-ACE+  and  Femstress. 

3.  The  MR  procedure  extracts  M  and  a  (partly  implicitly)  linearized  K  matrices  from  the  full 
FEM  solution,  at  the  applied  bias  Vo.  K  and  M  are  not  required  explicitly,  only  their 
product  with  arbitrary  vectors. 
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4.  The  MR  procedure  uses  Amoldi  Method  to  reduce  K  and  M  to  smaller  matrices  K'  and 
M'.  If  n  is  the  order  of  the  reduced  model,  i.e.  K'  is  (n  x  n)  and  M'  is  (n  x  n),  the  reduced 
ODE  system  matches  exactly  n  moments  of  the  full  problem. 

5.  Use  Matlab,  or  similar  matrix  solver,  to  solve  the  reduced  system  to  obtain  quickly  a 
frequency  characteristic  or  time-domain  (transient)  response. 

This  process  has  to  be  repeated  for  any  new  bias  point  (voltage  Vo). 

INPUT  for  the  MR  Procedure  (to  be  provided  by  user): 

-  name(s)  of  the  surface(s)  in  the  roM  model  where  the  voltage  is  applied; 

-  node  number  of  which  position/displacement  is  of  interest  as  a  result; 

-  defined  direction  (x,  y,  z)  of  the  node  displacement  of  interest; 

-  order  of  the  reduced  matrices. 

OUTPUT  from  the  MR  Procedure: 

-  the  reduced  matrices  in  the  input  format  (.in  file)  ready  for  Matlab  solution; 

-  sample  Matlab  scripts  (.m  files)  to  obtain  frequency  characteristic  or  transient  response. 


The  above  procedure  was  demonstrated  and  tested  on  the  example  of  a  MEMS  beam  actuated 
electrostatically  (Figure  145).  A  full  FEM  model  was  build  using  CFD-Micromesh  software 
from  a  layout.  CFD-ACE+  (Femstress)  was  used  as  a  coupled  electrostatic-stress  solver  for  the 
full  FEM  model. 


Figure  145.  A  coupled  electromechanical  system  used  as  an  example  for  the  model 
reduction  demonstration  with  the  Arnoldi  method. 
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Reduced  Model  vs.  Full  Model: 

-  full  FEM  model  order:  954  DOFs 

-  reduced  model  order:  10  DOFs 

-  CPU  time  acceleration  for  transient  simulations:  1  million  times  faster  per  time  step! 

Figure  146  shows  frequency-domain  characteristics  of  the  electrostatically  actuated  beam, 
calculated  with  different  methods:  full  FEM  non-linear  transient  simulation  (stars  in  the  plot), 
linearized  modal  analysis  in  CFD-ACE+  (full  line),  and  the  Amoldi  reduced  model  (circles). 

Frequency  Response  of  Beam  at  bias  =  100  v 


Figure  146.  Frequency  response  of  the  MEMS  beam  calculated  with  different  methods. 

As  it  can  be  seen  in  this  figure,  all  the  characteristics  match  very  well.  On  the  other  hand,  the  use 
of  the  reduced  model  in  a  transient  simulation  of  the  MEMS  device  (for  example,  with  Matlab) 
may  decrease  the  computation  time  by  a  factor  of  1  million  per  one  time  step!  This  is  very 
promising  for  system-level  simulations  of  large  MEMS  circuits. 

6  CONCLUSIONS  AND  RECOMMENDATIONS 

6.1  Conclusions  from  the  Project 

This  report  has  presented  the  current  status  and  results  of  the  DARPA  BAA  97-17  project  aimed 
at  developing,  demonstration,  and  validation  of  new  concepts  of  reduced/compact  models  of 
electro-mechanical,  fluidic,  and  thermal  phenomena  in  MEMS,  as  well  as  at  development  of  new 
Software  System  (in  CFD-ACE-i-)  for  Automated  Generation  of  Reduced  or  Compact  Models  of 
Microdevices  from  High  Fidelity  3-D  Simulations.  All  of  the  objectives  of  the  project  were  met 
successfully.  Many  new  reduced  model  concept  have  been  developed  and  validated,  both 
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through  comparison  with  high-fidelity  simulations  and  experimentally.  New  model  generation 
procedures  have  been  proposed,  implemented,  and  tested.  New  reduced  models  have  been 
demonstrated  in  output  formats  compatible  with  system-level  simulators  like  SPICE  or 
Saber/MAST.  The  air  damping  compact  models  resulting  from  this  project  have  been  included 
into  NODAS  library  from  CMU  for  the  MEMS  system-level  design. 

Several  new  interfaces  and  import/export  procedures  for  CAD/EDA  standards  have  been  added 
to  CFD-ACE-t-  software  package,  to  enable  easy  and  automatic  building  of  3-D  models  from 
EDA  layouts  and  process  data.  The  CFD-ACE-i-  solver  itself  has  been  extended  and  improved  for 
reliable  high-fidelity  simulations  of  coupled  physical  phenomena  specific  for  microdevices.  New 
algorithms  and  tools  have  been  implemented  in  CFD-ACE-t-  to  enable  the  use  of  reduced  or 
mixed  dimensionality  in  ID,  2D,  or  3-D  simulations  of  MEMS  devices  and  systems. 

The  major  achievements  of  the  program  have  already  been  outlined  before,  and  are  presented 
below  in  brief: 


New  Reduced/Compact  Models  of  Microdevices 

•  New  nonlinear  compact  models  of  squeeze  film  damping  in  MEMS,  validated  with  CFD- 
ACE-f  3-D  simulations,  and  implemented  in  SPICE  and  SABER  formats,  and 
implemented  in  NODAS  MEMS-CAD  system  at  CMU. 

•  The  new  compact  squeeze  film  model  was  successfully  compared  with  experimental  data 
from  a  CMOS-MEMS  bandpass  filter.  The  NODAS  simulation  results  with  the  new 
damping  model  are  in  excellent  agreement  with  the  measurements. 

•  New  compact  model  of  lateral  (shear)  damping  in  MEMS,  validated  with  CFD-ACE-i-. 

•  Developed  mixed-dimensional  and  reduced  fluidic  models  of  microchannels, 
microvalves,  micropumps,  droplet  generators,  synthetic  jets,  and  piezoelectric 
micropump. 

•  Novel  concept  of  compact  models  for  synthetic  jets,  using  a  polyhedral  control  volume 
capability  of  CFD-ACE-f.  The  model  was  validated  against  3-D  high-fidelity  simulation 
data  obtained  for  a  range  of  parameters. 

•  New  equivalent-circuit  models,  in  SPICE  and  SABER  formats,  for  arbitrary-shaped 
fluidic  microchannels,  and  procedures  of  their  generation. 

•  New  Matrix  Reduction  Techniques  for  structure  dynamics:  Lanczos  method  and  WYD 
Ritz  algorithm,  developed  at  UF.  Tests  showed  up  to  10,000x  acceleration  in  comparison 
to  full  3-D  FEM  computation  time. 

Interfaces  And  Import/Export  Procedures  for  CAD/EDA  Standards.  To  Enable  Easy  And 

Automatic  Building  of  3-D  Models  from  EDA  Layouts  And  Process  Data 

•  Enhancement  of  CFDRC  software  to  process  CAD  data  for  MEMS  design.  We  are  now 
able  to  read  input  files  in  the  formats:  DXF,  IGES,  PATRAN,  CIF,  GDSII,  and  several 
others. 

•  Two  new  data  formats,  CIF  and  GDSII,  have  been  implemented  in  CFD-GEOM, 
enabling  input  and  conversion  from  external  files  coming  from  EDA  tools. 

•  Implementation  of  VRML  (Virtual  Reality  Modeling  Language)  into  the  computational 
environment  of  CFDRC. 
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•  A  new  software,  CFD-Micromesh,  for  automatic  3-D  solid  model  generation  and 
meshing  from  MEMS  layouts  imported  directly  from  CAD/EDA  systems.  The  CIF  and 
GDSII  layout  files  are  imported  easily  into  CFD-Micromesh,  allowing  coupling  the 
software  with  several  commercial  IC  design  tools.  The  new  tool  is  very  fast  and  equipped 
with  user-friendly  graphical  user  interface.  By  clicking  only  a  single  button  in  CFD- 
Micromesh,  a  three-dimensional  solid  model  is  generated  fully  automatically  from  layout. 
Building  the  3-D  model  as  well  as  automatic  generation  of  a  full  3-D  unstructured 
computational  mesh,  again  only  by  clicking  a  single  button  in  CFD-Micromesh,  is 
performed  typically  within  several  minutes  on  a  current  PC  workstation.  The 
automatically  generated  three-dimensional  device  model,  with  the  unstructured  3-D 
computational  mesh  in  DTF  format,  is  imported  directly  into  CFDRC's  simulator  CFD- 
ACE+. 

•  Developed  interfaces  between  CFDRC  software  and  Cadence  design  tools. 

•  A  new  data  format  (called  GBV  -  Geometry,  Boundary  and  Volume  conditions)  for 
interchange  of  simulation  data  between  Cadence  Design  Framework  II  schematic  view 
and  CFDRC  ACE+  simulator  has  been  developed  at  CMU. 

New  Features  in  CFD-ACE-h  for  High-Fidelity  Simulations  Specific  for  MEMS 

•  Second  order  boundary  conditions  for  wall  treatment. 

•  Second  order  time  treatment,  through  modified  Crank-Nicolson  scheme. 

•  Automatic  grid  remeshing  for  arbitrary  moving  boundaries/bodies  in  6  degrees  of 
freedom  (DOFs).  A  trans-finite  interpolation  (TFI)  technique  with  mesh  smoothing  was 
developed  and  implemented,  which  allows  for  effective  simulation  of  dynamic  3-D 
motion  of  mirrors,  membranes,  gyroscopes,  and  other  MEMS  structures. 

•  New  capability  in  CFD-ACE+  for  extraction  lumped  parameters  (normal/tangential  force, 
torque,  total  pressure  and/or  flow)  for  compact  model  generation. 

•  "Property  Manager"  -  new  interactive  framework  software  for  MEMS  material  property 
database. 

Procedures  And  User-Friendly  Interfaces  For  Extraction  Of  Compact-Model  Parameters  and 

Characteristics  From  High-Fidelity  Simulations 

•  Capability  to  perform  parametric  runs  in  CFD-ACE-t-  through  graphical  user  interface,  for 
mechanical,  electrostatic,  and  fluidic  simulations,  both  steady-state  and  transient. 

•  A  dedicated  computer  program  Squeeze,  written  in  portable  C++  language,  for  automatic 
generation  of  input  list  in  SABER  or  SPICE  format,  describing  equivalent-circuit  model 
of  squeeze-film  damping  between  moving  plates. 

•  Procedures  for  automatic  generation  of  compact  model  (extraction  of  lumped 
capacitances)  of  a  comb-drive  resonator,  using  results  of  high-fidelity  3-D  electrostatic 
simulations  with  FastBEM  from  CFDRC. 

•  Procedures  for  automatic  generation  of  equivalent-circuit  models  of  microfluidic  devices 
from  parametric  runs  of  CFD-ACE+;  examples  of  SPICE  and  SABER  compact  models 
of  Tesla  valve. 

•  Procedures  of  FEM  model  reduction  based  on  Amoldi  method  for  electro-mechanical 
microdevices,  developed  at  Prof.  Jacob  White’s  group  at  MIT,  have  been  implemented 
and  installed  at  CFDRC  and  demonstrated  and  tested  in  coupling  with  CFD-ACE+  code. 


136 


Libraries  Of  Reduced  Models  for  MEMS  CAD  System-Level  Tools 


•  The  new  compact  models  of  squeeze  film  damping  have  been  incorporated  into  the 
NODAS  MEMS-CAD  system  library  from  CMU. 

•  The  new  compact  models  of  lateral  shear  damping  have  also  been  inserted  into  the 
NODAS  MEMS-CAD  system  library. 

•  The  frequency  dependent  transfer  function  of  lateral  damping  has  been  implemented  in 
Verilog-A  as  a  Laplace  transform  and  used  in  NODAS. 

Validation/Demonstration  of  the  New  CFD-ACE-(-  Features  and  the  New  Reduced  Models 

•  At  Georgia  Tech  (GT); 

Fabricated  micromachined  synthetic  jets  with  integrated  MEMS  modulators. 

-  Several  methods  have  been  used  to  test  the  synthetic  jets  and  modulator  arrays:  visual 
deflection  test,  flow  modulation  test,  Schlieren  visualization,  pitot  tube,  x-wire 
anemometer,  and  particle  image  velocimetry  (PIV). 

Characteristics  of  performance  of  arrays  of  synthetic  jets  in  local  active  control  of 
microsystem  temperature  were  measured  with  infra-red  (IR)  camera. 

•  The  PrV  flow  data  from  GT  were  compared  successfully  with  2D  and  3-D  numerical 
simulations  from  CFD-ACE+. 

•  The  temperature-system  data  measured  at  GT  with  ER  camera  were  compared 
suecessfully  with  numerical  simulations  at  CFDRC,  including  the  use  of  reduced  models 
in  jet  arrays. 

•  At  University  of  Florida  (UF): 

Validation  of  the  implementation  of  Lanczos  algorithm  through  comparison  with  a 
MEMS  structure  from  MIT. 

•  At  Carnegie  Mellon  University  (CMU): 

-  Experimentally  verified  the  squeeze  film  model  using  a  three  resonator  CMOS 
micromachined  mechanical  filter  to  show  that  the  improved  squeeze  film  model 
reduced  error  in  quality  factor  (Q)  computation  from  20%  to  2%. 

Experimentally  verified  the  lateral  damping  model  using  folded  flexure  MUMPS 
resonators.  Inclusion  of  finite  size  effects  reduced  errors  from  20%  to  8%. 

•  At  CFDRC,  both  the  high-fidelity  3-D  CFD  simulations  of  Tesla  valve  and  its  compact 
models  in  SPICE  and  Saber  were  successfully  compared  with  characteristics  measured  at 
University  of  Washington. 

The  mixed-dimensionality,  multi-energy  domain  capability  in  CFD-ACE-l-  is  the  first  ever 
computational  software  applicable  to  both  micro-scale,  system-scale,  and  mixed  micro/system- 
scale  simulations.  The  macro-scale  (package,  board)  allows  the  user  to  verify  the  performance  of 
the  microdevice  (sensor,  transducer,  actuator,  power  device,  optoelectronic  array)  in  an  external 
system  or  environment. 
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6.2 


Commercial  Effects  of  the  Project 


6.2.1  "Success  Stories" 

•  CFDRC  has  been  awarded  a  new  SBIR  contract  in  2000,  for  U.S.  Air  Force 
(AFRL/VSSE),  to  develop  new  “High  Speed  Inertial  MEMS  Sensors  with  Field  Emitter  Array 
Readout”  for  application  in  novel  gyroscopes  for  satellite  and  vehicle  positioning,  etc. 

•  CFDRC's  High-Fidelity  and  Compact  Models  of  Ink-Jets  used  by  Kodak. 

•  Reduced  Models  of  Synthetic  Jets  used  for  Virtual  Flight  Control  Simulation. 

•  Air-Damping  Compact  Models  under  development  for  Lucent  Technologies,  for  analysis 
and  design  of  MEMS-based  micromirrors  for  fast  optical  switching. 

•  Reduced  Models  of  Synthetic  Jets  will  be  used  within  the  DARPA  HERETIC  Program 
(Heat  Removal  by  Thermo-Integrated  Circuits),  where  CFDRC  will  support  modeling  and 
design  of  Arrays  of  Synthetic  Microjets  for  active  cooling  of: 

VCSEL  Arrays  (Linear  and  2D)  Georgia  Tech  +  Honeywell 

Microprocessor  Chips  -  Georgia  Tech  +  Intel. 

•  Air-Damping  Compact  Models  used  by  Delphi  Automotive  Systems. 

The  last  point  is  presented  below  in  more  detail. 

Air-Damping  Compact  Models  used  by  Delphi  Automotive  Systems 

Compact  models  of  air  damping  for  MEMS  devices,  developed  by  CfDRC  under  DARPA 

funding  within  the  Composite  CAD  Program,  appeared  very  helpful  in  the  analysis  process  of 

accelerometers  for  air-bag  deployment  systems,  being  developed  by  Delphi  Automotive 

Systems.  The  squeeze-film  compact  models  allow  to  model  the  system  using  Saber  simulator. 


•  A  compact  model  of  squeeze-film  damping  effects  acting  on  accelerometer  moving  plate 
was  developed  in  Saber  (MAST)  format,  using  three-dimensional  CFD-ACE-l- 
simulations  in  frequency  domain  for  parameter  extraction. 


For  the  model  validation,  a  wider  range  of  operating  conditions  (for  instance,  larger 
amplitudes,  etc.)  has  been  applied,  and  the  results  obtained  in  Saber  using  the  compact 
model  were  compared  to  results  of  three-dimensional  simulations  of  CFD-ACE+. 
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3-D  Simulation  -  CFD-ACE+ 


System-level  simulation  -  SABER 
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•  The  reduced  model  of  squeeze-film  damping  appeared  to  be  very  helpful,  both  for  Delphi 
and  CFDRC,  in  understanding  some  physical  phenomena  involved,  and  in  the  analysis  of 
the  accelerometer  structure  behavior. 

Success  of  CFDRC  in  Air  Damping  Modeling  and  Analysis  in  MEMS 

As  a  result  of  CFDRC  activities  in  this  project,  and  resulting  publications  and  presentations, 
CFDRC  was  approached  by  numerous  commercial  and  academic  institutions,  both  from  USA 
and  international,  which  expressed  an  interest  in  methods  and  tools  for  modeling  the  air  damping 
phenomena  in  MEMS.  It  appears  that  thanks  to  this  program  CFDRC  has  become  the  deHnite 
world  leader  in  capabilities  of  viscous  damping  modeling  in  MEMS,  both  in  terms  of  high- 
fidelity  flow  modeling  as  well  as  in  the  generation  and  validation  of  reduced/compact  models. 

Institutions,  which  approached  CFDRC  with  interest  in  Air  Damping  modeling  and  analysis  in 
MEMS,  including  Compact  Models,  are  as  follows: 


Name 

Institution 

u.s. 

1. 

Arjun  Selvakumar 

Input/Output,  Inc.,  USA 

2. 

Bill  Grande 

Kodak  Research  Labs,  Rochester,  NY,  USA 

3. 

Brady  Davies 

Kistler  Instrument  Corporation,  USA 

4. 

Calvin  Adkins 

Harris  Corporation,  Melbourne,  FL,  USA 

5. 

Clarence  Chui 

Mdigm  Display,  San  Francisco,  CA,  USA 

6. 

Dan  Koch 

Motorola,  Semiconductor  Products  Sector,  USA 

7. 

David  Rich 

Delphi  Automotive  Systems,  Sensor  Design,  Kokomo,  IN,  USA 

8. 

Dick  Nelson 

MCC,  Austin,  TX,  USA 

9. 

Edward  Chan 

Bell  Labs,  Lucent  Technologies,  NJ,  USA 

10. 

Edward  S.  Kolesar 

Texas  Christian  University,  Dept  of  Eng.,  Fort  Worth,  TX,  USA 
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11. 

Heinz  Busta 

Samoff  Corp.,  Princeton,  NJ,  USA 

12. 

Jason  Tauscher 

Silicon  Designs,  Inc.,  USA 

13. 

Joe  Seeger 

UC  Berkeley,  EECS  Department,  CA,  USA 

14. 

Kamil  Ekinci 

Condensed  Matter  Physics,  CalTech,  Pasadena,  CA,  USA 

15. 

Kophu  Chiang 

Tellium,  Inc.,  Oceanport,  NJ,  USA 

16. 

Richard  Johansen 

Rensselaer  Polytechnic  Institute,  Troy,  NY,  USA 

17. 

Robert  Conant 

UC  Berkeley,  CA,  USA 

18. 

Robert  S.  Okojie 

Ford  Microelectronics,  Colorado  Springs,  CO,  USA 

19. 

Sitaraman  V.  Iyer 

ECE  Dept,  Carnegie  Mellon  University,  Pittsburgh,  PA,  USA 

International 

20. 

Cho  NamKyu 

Korea  Electronics  Technology  Institute,  Korea 

21. 

Eskild  R.  Westby 

SensoNorasa,  Norway 

22. 

George  Raicevich 

National  Acoustic  Laboratories,  Chatswood,  NSW,  Australia 

23. 

Karsten  Funk 

Robert  Bosch  GmbH,  Germany 

24. 

Kazuyuki  Minami 

Dept  of  Mechatronics  and  Precision  Eng.,  Tohoku  Univ.,  Sendai, 
Japan 

25. 

Matthias  Maute 

Sensortechnologiezentrum,  Robert  Bosch  GmbH,  Reutlingen, 
Germany 

26. 

Patrick  Aldebert 

SUPELEC,  Service  des  Mesures,  Gif/Yvette,  France 

27. 

Peter  Frere 

LucasVarity  Automotive  Tech.  Centre  (part  of  TRW  Automotive), 
UK 

28. 

Rose  Zhang 

Photonics  Mat.  &  Proc.,  National  Optics  Institute,  Que,  Canada 

29. 

Roumiana  Paneva 

X-FAB  GmbH,  Erfurt,  Germany 

30. 

Roy  Knechtel 

X-FAB  GmbH,  Erfurt,  Germany 

31. 

Sangwoo  Lee 

Seoul  National  University,  EE  Dept,  Korea 

32. 

Shawn  Taylor 

Virtual  Science  Limited,  West  Sussex,  UK 

33. 

Tiansheng  Zhou 

Alberta  Microelectronic  Corporation,  Edmonton,  Alberta,  Canada 

34. 

Youngho  Kim 

Hanyang  University,  Seoul,  Korea 

CFD-Micromesh:  Fast  3-D  Model  Builder  from  Layouts 

The  new  CFD-Micromesh  software,  developed  during  this  project,  for  automatic  3-D  solid 
model  generation  and  meshing  from  electronic  and  MEMS  layouts,  has  a  very  big  and 
promising  commercial  potential.  The  popular  layout  description  formats,  CIF  and  GDSII,  are 
imported  easily  into  CFD-Micromesh,  allowing  coupling  the  software  with  several  commercial 
EDA  design  tools.  The  model  building  in  CFD-Micromesh  is  very  fast  and  controlled  by  a  user- 
friendly  graphical  user  interface  (GUI).  The  automatically  generated  three-dimensional  device 
models,  using  the  unstructured  3-D  computational  mesh  in  DTF  format,  can  be  imported  directly 
into  CFDRC's  simulator  CFD-ACE-i-  for  multi-physics  simulations. 

CFD-Micromesh  was  recently  extended  with  a  capability  to  read  a  fabrication  process  (layer 
structures)  description  in  a  standard  format  called  SIPPs  (Standard  Interconnect  Performance 
Parameters).  The  format  is  being  developed  by  several  big  semiconductor  companies,  like  Si2, 
AMD,  Motorola,  IBM,  Cadence,  Avanti,  Mentor  Graphics,  and  others,  with  participation  of 
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CFDRC  (more  info:  http://www.si2.org/sipps/  ).  SIPPs  layer  data  are  typically  imported  in 
connection  with  a  layout  in  GDSn  format,  which  makes  together  a  complete  3-D  description  of  a 
design.  This  capability  makes  CFD-Micromesh  a  very  attractive  modeling  tool  also  for  VLSI 
electronics  industry,  in  particular  for  analysis  of  interconnects  and  PCBs.  CFDRC  has  already 
started  an  aggressive  entering  into  this  market,  and  we  are  planning  to  expand  it  in  near  future. 

As  a  part  of  the  marketing  campaign,  a  free  evaluation  version  of  CFD-Micromesh  is  offered 
for  download  (for  a  limited  time)  from  the  web  site  http://www.cfdrc.com/~micromesh. 

This  version  of  CFD-Micromesh  imports  CIF  and  GDSII  layouts,  reads  the  SIPPs  process 
description  format,  and  includes  direct  interface  to  Virtuoso^™^  layout  editor  from  Cadence.  The 
evaluation  version  builds  automatically  3-D  solid  models  from  layout  and  process  data,  and 
allows  the  user  to  view  them  interactively. 

Finally,  CFDRC  has  negotiated  an  agreement  with  Cadence  Design  Systems,  Inc.  to  develop 
interfaces  and  integration  of  CFDRC's  MEMS  simulation  software  with  Cadence  CAD  tools. 
CFD-Micromesh  will  be  a  part  of  Cadence  CAD  Environment,  used  as  a  “3-D  Project  Viewer” 
from  Virtuoso  layout  editor  from  Cadence  (see  Figure  5). 

6.2.2  Technology  Transfer 
Technology  Transfer  to  Industry 

Additionally  to  the  industrial  contacts  and  projects  mentioned  in  the  previous  section,  as  a  result 
of  this  project  CFDRC  provided  software,  training,  and  support  oriented  for  MEMS  design  to  the 
following  companies: 

•  Aclara 

•  Procter  and  Gamble 

•  Kodak 

•  Lucent  Technologies 

•  Hewlett-Packard 

Technology  Transfer  to  Academia 

The  academic  parties  of  this  project  have  used  CFD-ACE+  software  as  a  design,  research,  and 
educational  tool.  CFDRC  has  provided  software,  training,  and  support  in  exchange  for 
constructive  feedback,  technology  qualification,  and  recommendations.  The  following  academic 
institutions  have  been  involved  in  those  of  activities  within  this  project: 

•  Carnegie  Mellon  University 

•  University  of  Florida 

•  Stanford  University 

•  Cornell  University 

•  Georgia  Tech 

•  University  of  California  Los  Angeles  (UCLA) 

•  University  of  California  Berkeley 

•  Drexel  University 
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•  University  of  California,  Irvine 

•  NC  State  University 

•  Pennsyl.  State  University 

•  MIT 

•  University  of  South  AL 

•  University  of  Central  FL 

•  University  of  Sydney 

•  University  of  Washington 

•  University  of  Neuchatel 

•  University  of  Minnesota 

Talks  /  Presentations  (not  listed  in  the  published  References) 

The  progress  of  the  project  was  presented  at  each  DARPA  Composite  CAD  PI  Meeting,  between 

October  1997  and  August  2000.  In  addition,  the  following  presentations  related  to  this  project 

took  place: 

•  During  the  MSM'98  Conference  (April  6-8,  1998,  Santa  Clara,  CA),  three  presentations 
were  made: 

A.J.  Przekwas  and  A.  Krishnan,  “ACE+MEMS:  An  Integrated  Multi-Disciplinary 
CAD  System  for  Micro-Electro-Mechanical  Systems  (MEMS)” 

H.Q.  Yang  and  A.J.  Przekwas,  “Computational  Modeling  of  Microfluidic  Devices 
with  Free  Surface  Liquid  Handling” 

M.M.  Athavale,  H.Y.  Li,  A.J.  Przekwas,  B.  Piekarski,  R.  Zebo,  and  S.  Tenney, 
“Modeling  3-D  Fluid  Flow  for  a  MEMS  Laminar  Proportional  Amplifier”  (in 
collaboration  with  Army  Research  Laboratory,  Adelphi,  MD) 

•  At  the  American  Institute  of  Chemical  Engineering  (AIChE)  Conference  on 
Microreactors  in  New  Orleans,  Louisiana,  March  8-12,  1998,  a  CFDRC  paper  entitled, 
“Integrated  Computational  Tools  for  Interdisciplinary  Design  of  Microreactors”,  by  M.G. 
Giridharan,  A.  Krishnan,  S.  Krishnamoorthy,  and  A.J.  Przekwas,  was  presented. 

•  Dr.  A.J.  Przekwas  was  invited  to  the  European  Workshop  on  MEMS  Modeling 
(organized  by  Prof.  J.  Korvink  of  the  Univ.  of  Freiburg,  Germany,  in  fall  1998)  to  present 
ACE-i-MEMS  R&D  projects  and  software. 

•  Presentation  "Generation  of  Reduced  Models  for  MEMS"  (by  Andrzej  Przekwas, 
CFDRC)  at  the  Workshop  on  Interdisciplinary  Design  and  Simulation  Tools  for 
Microsystems,  in  Munich,  Germany,  June,  1998. 

•  Presentation  "Reduced  Model  Generation  for  Microfluidic  Systems"  (by  Andreas  Klein 
and  Andrzej  Przekwas,  CFDRC)  at  the  2nd  Workshop  on  Interdisciplinary  Design  and 
Simulation  Tools  for  Micro-  and  Biomedical  Fluidic  Applications  in  Berlin,  Germany, 
June  17,  1999. 

•  “Multi-Physics  Simulation  Tools  For  Design  And  Optimization  Of  CBW  Sensors”, 
presented  by  Andrzej  J.  Przekwas  and  Vinod  B.  Makhijani,  at  DARPA  SIMBAD 
Meeting,  April  10,  2000,  Washington  DC 
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6.3  Recommendations  and  Plans  for  Further  Research 


During  the  course  of  this  project  several  new  capabilities  have  been  implemented  in  CFD-ACE+ 
facilitating  automated  creation  of  3-D  models  and  automated  generation  of  compact  models  for 
MEMS  devices  for  system  level  simulations.  In  the  last  phase  of  the  project,  focused  on 
demonstrations  and  practical  applications,  we  received  a  substantial  feedback  and  assessment  of 
the  ACE+  code  in  the  area  of  multiphysics  modeling,  model  setup,  interfaces  with  EDA  tools, 
compact  model  generation,  and  engineering  data  extraction.  ACE+  proved  robust  in  generation 
of  reduced  models  for  linear  problems  e.g.  electrostatics,  deformation,  and  unique  for  nonlinear 
air  dumping.  Several  customers  requested  macromodel  extraction  capability  for  coupled  domain 
physics  e.g.  electrostatics  and  deformation,  thermal  and  mechanical  disciplines.  We  have 
developed,  demonstrated,  and  delivered  this  capability  in  the  last  year.  Despite  several  efforts  of 
our  team  and  others  (MU,  CMU,  ...)  no  robust  general  method  could  be  established  for 
nonlinear  or  strongly  coupled  problems  e.g.  fluidics,  mixing,  chemistry. 

For  nonlinear  problems  we  recommend,  and  help,  our  customers  to  use  simplified  mathematical 
formulation  of  general  conservation  equations  with  semi-empirical  parameters  derived  from  3-D 
simulations  e.g.  synthetic  jet  models,  orifices,  valves,  . .  .We  have  also  developed  a  new  modeling 
concept  “Filament  Model”  to  facilitate  high  (system)  level  modeling  in  ACE+.  This  technique  is 
being  further  developed  in  the  parallel  DARPA  project  on  “VLSI  CAD  for  Microfluidics”. 
CFDRC  is  fully  committed  to  further  development  of  system-level  simulations  and  model 
extraction  for  microsystems.  The  demand  for  this  capability  will  be  steady  growing  as  the  basic 
MEMS  and  microfluidic  devices  are  better  established.  We  will  further  evolve,  improve,  and 
commercialize  the  model  extraction  and  system-level  modeling  with  existing  resources  from 
CFDRC,  government,  and  industry. 

Based  on  the  experience  gained  in  this  project  we  have  identified  several  areas  for  further 
research,  model  development,  and  CFD-ACE+  software  improvement.  They  can  be  grouped  in 
following  categories: 

A.  Interfaces  with  EDA  Tools 

In  the  present  project,  we  have  developed  new  tool,  CFD-Micromesh,  to  interface 
microelectronic  and  mechanical  CAD  domains.  The  code  is  rapidly  gaining  acceptance  and 
interest  from  broad  range  of  applications.  It  is  one  of  the  biggest  achievements  of  this  project.  It 
has  however  several  limitations  for  applications  in  microfluidic  MEMS  and  for  BioMEMS.  The 
most  urgent  CFD-ACE-i-  needs  for  system-level  design  are: 

•  semiautomatic  block  structured  grid  generation, 

•  mixed  structured-unstructured  grids, 

•  adaptive  octree  grids, 

•  two-way  link  with  CFD-GEOM, 

•  parametric  modeling  in  Micromesh, 

•  Python  scripting  model  setup  using  Micromesh  functions, 

•  automated  construction  of  filaments  in  Micromesh, 

•  tighter  link  with  process  fabrication  and  layout  tools, 

•  boolean  operations  on  masks  for  multilayer  microstructures. 
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B.  Automation  and  Parameterization  of  Model  Setup 

At  present,  CFD-ACE+  is  suitable  for  single  runs  of  a  specific  design  or  pre-arranged  multi-runs 
for  single  variation  in  e.g.  boundary  conditions.  There  is  a  strong  need  for  parametric  modeling, 
scripting,  optimization,  and  inverse  problem  solutions.  A  first  step  has  been  recently  made  at 
CFDRC  by  developing  a  Python  scripting  capability  for  CFD-GEOM.  We  hope  to  expand  this 
capability  to  all  other  modules  of  ACE+,  like  GUI,  DTE,  VIEW,  and  Micromesh.  Advanced 
academic  and  industrial  users  should  be  able  to  link  the  functions  and  modules  in  ACE-t- 
modeling  dataflow.  Coupling  of  ACE-t-  with  an  optimization  module  would  answer  requests  from 
several  customers  who  have  requested  device  optimization  capability. 

C.  Mixed-Dimensionality  Modeling  for  Micro-Bio-Fluidics 

In  several  circumstances,  it  will  be  very  difficult  to  develop  compact  or  circuit  models  for  key 
devices,  like  bioreactor,  biodetector,  etc.  We  would  like  to  develop  a  mixed-level  and  mixed- 
dimensionality  modeling  capability.  We  have  demonstrated  it  on  the  coupling  of  ACE+  with 
SABER.  There  is  a  big  academic  and  commercial  interest  to  develop  direct  interface  between 
ACE+  and  SPICE.  We  also  plan  to  expand  the  Filament  modeling  concept  for  mixed- 
dimensionality  modeling.  Further  improvements  in  model  setup,  meshing,  and  data  post¬ 
processing  are  needed. 

D.  Extension  of  Compact  Model  Generation  for  New  Microdevices 

New  DARPA  programs  have  been  initiated  in  BioChips,  Power  MEMS,  and  other  new  areas  are 
anticipated,  e.g.  MOEMS,  RF-MEMS,  Bio-Photonics,  bio-medical,  etc.  Existing  compact  model 
extraction  capabilities  in  ACE-i-  have  been  shown  for  mechanical,  thermal,  and  some  fluidic 
devices.  There  is  a  strong  demand  (from  other  DARPA  project  teams  in  particular)  for  compact 
models  for  other  disciplines  such  as  electrokinetics,  electrochemistry,  biochemistry,  photonic 
devices,  optoelectronic  devices,  VLSI  interconnects,  fluidic  cooling  circuits,  etc. 

E.  Expansion  of  Model  Libraries 

Most  of  the  high-fidelity  models  constructed  in  this  project  are  nonparametric  models  saved  as 
DTE  files.  There  is  an  absolute  demand  for  parametric  models.  We  have  started  to  develop 
Python  scripts  for  MEMS.  There  is  a  great  potential,  but  a  lot  remains  to  be  done  to  implement 
Python  in  all  other  tools. 

We  would  like  to  offer  Python  scripts  over  the  Internet,  and  allow  other  users  (e.g.  students)  to 
upload  their  models  and  share  them.  We  will  also  explore  the  idea  of  e-business  with  others. 
Closer  collaboration  with  the  MEMS  Exchange,  which  has  already  started,  would  be  beneficial. 

F.  Strong  Interface  between  ACE+  and  System  Level  CAD  (SPICE) 

Most  of  the  system-level  simulations  in  the  Composite  CAD  Program  were  tested  and 
demonstrated  with  the  Saber  software  from  Analogy.  The  Analogy  company  has  been  recently 
acquired  by  Avant!,  and  the  future  of  Saber  and  its  market  position  is  not  sure.  Gn  the  other 
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hand,  behavioral  models  implemented  in  Saber  had  advantage  over  the  SPICE  models  in  that 
they  offered  more  direct  representation  and  simulation  of  non-electrical  phenomena  (fluidic, 
mechanical,  thermal,  etc),  readily  available  to  designers  in  hydraulic  or  mechanical 
displacement/force  units,  without  additional  translation  to  electrical  equivalents.  The  MAST 
language  in  Saber  has  also  offered  a  very  flexible  way  of  adding  to  simulations  custom-designed 
models  generated  by  users.  This  is  why  CFDRC  has  developed  new  software  tools  in  CORBA 
to  enable  direct  coupling  between  Saber  and  ACE+  simulations. 

Nevertheless,  SPICE  remains  the  most  popular  circuit  simulator  in  the  world,  and  it  is  available 
in  several  free  versions  for  many  users.  During  this  project,  we  have  demonstrated  generation  of 
a  compact  model  of  an  arbitrary  microfluidic  channel,  in  the  form  of  nonlinear  equivalent  circuit. 
It  is  very  attractive  and  useful  to  a  wide  community  of  SPICE  users  who  are  very  often  also 
designers  of  microfluidic  MEMS  and  microsystems  in  the  domain  of  bio-medical  and  chemical 
engineering,  integrated  very  often  with  control  and  drive  electronics.  Because  of  that,  and  also 
due  to  not  sure  future  of  Saber,  CFDRC  will  pursue  the  idea  of  integrating  ACE+  with  SPICE 
circuit  modeling  capabilities.  In  this  endeavor  we  plan  to  collaborate  with  the  group  of  Michael 
Shur  at  RPI  (developers  of  AIM-SPICE),  and  the  UC  Berkeley  CAD  Group  -  original  creators  of 
SPICE  and  CIDER,  as  well  as  UC  Berkeley  BSAC  (Kris  Pister)  Group  which  developed 
SUGAR,  a  system-level  simulation  program  for  MEMS  devices. 

G.  Software  Verification  on  Design  of  Practical  Microsystems 

CFDRC  is  involved  in  MEMS  related  projects  with  several  academic,  commercial,  and 
government  partners.  A  lot  of  effort  is  provided  to  support  our  partners.  We  need  to  improve  and 
adapt  CFD-ACE+  environment  to  allow  more  flexible  access  to  the  modules  of  the  code  and  to 
the  dataflow. 

More  partnership  projects  with  US  industry  are  needed,  and  planned,  to  adapt  the  code  for 
practical  design  environments.  Much  effort  is  needed  to  simplify  and  automate  the  design 
process  for  non-experts.  Several  of  our  key  customers  indicated  that  there  is  no  available 
graduates  with  multiphysics  expertise  needed  to  operate  the  high-fidelity  software  such  as 
CFD-ACE+.  This  indicates  that  even  larger  effort  is  required  from  CjTORC  to  make  our 
multiphysics  software  more  accessible  for  new  users  and  to  allow  them  a  faster  learning  curve. 
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